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ABSTRACT 


Design  requirements  for  a  small  satellite  (NPSAT-1)  Attitude  Determination  and 
Control  Subsystem  (ADCS)  is  a  three-axis  stabilized  spacecraft  which  requires  a  control 
attitude  of  +/-  1.0  degrees  and  knowledge  attitude  of  +/-  0.1  degree.  Several  design 
aspects  are  considered  in  development  of  attitude  control  systems  for  a  small  satellite, 
such  as:  spacecraft  dynamics,  space  environment,  disturbance  torques,  orbit  type,  and 
spacecraft  complexity.  The  ideal  spacecraft's  attitude  sensor  is  a  rate  gyroscope,  which 
provides  rate  infonnation  to  the  attitude  control  system.  In  the  case  of  NPSAT-1,  due  to 
budget  constraints  alternative  sensors  will  be  utilized,  such  as:  a  three-axis 
magnetometer,  earth  sensors,  and  a  Global  Positioning  System  (GPS).  A  small  satellite 
designed  to  have  a  three-axis  stabilized,  biased  momentum  system,  must  have  a  robust 
control  system,  and  requires  a  momentum  wheel  to  provide  stiffness  to  maintain  attitude, 
and  magnetic  torque  rods  on  each  axis.  The  current  design  of  NPSAT-1  uses  all  of  these 
sensors  to  provide  rate  information  for  damping  and  stability  to  the  control  system  that 
requires  a  complicated  attitude  control  design.  The  purpose  of  this  attitude  control  design 
simulation  is  to  investigate  and  propose  a  control  law  utilizing  a  single  pitch  momentum 
wheel  and  three  magnetic  torque  rods.  A  further  proposal  is  to  utilize  a  constant  speed 
momentum  wheel  to  avoid  momentum  damping  and  over  speed,  replace  the  pitch  control 
with  magnetic  torquers,  and  develop  a  Kalman  filter  estimator  to  provide  all  the  required 
angular  rates. 
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DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed  in  this  simulation  may 
not  have  been  exercised  for  all  cases  of  interest.  While  every  effort  has  been  made, 
within  the  time  available,  to  ensure  that  the  programs  are  free  of  computational  and  logic 
errors,  they  cannot  be  considered  completely  validated.  Any  application  of  these 
programs  without  additional  verification  is  at  the  risk  of  the  user. 
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EXECUTIVE  SUMMARY 


Design  requirements  for  a  small  satellite  (NPSAT-1)  Attitude  Determination  and 
Control  Subsystem  (ADCS)  is  a  three-axis  stabilized  spacecraft  that  requires  a  control 
attitude  of  +/-  1.0  degrees  and  knowledge  attitude  of  +/-  0.1  degree.  Several  design 
aspects  are  considered  in  development  of  attitude  control  systems  for  a  small  satellite: 
spacecraft  dynamics,  space  environment,  disturbance  torques,  orbit  type,  and  spacecraft 
complexity.  The  ideal  spacecraft's  attitude  sensor  is  a  rate  gyroscope,  which  provides 
rate  information  to  the  attitude  control  system.  In  the  case  of  NPSAT-1,  due  to  budget 
constraints  alternative  sensors  will  be  utilized:  a  three-axis  magnetometer,  earth  sensors, 
and  a  Global  Positioning  System  (GPS).  A  small  satellite  designed  to  have  a  three-axis 
stabilized,  biased  momentum  system  must  have  a  robust  control  system  and  requires  a 
momentum  wheel  to  provide  stiffness  to  maintain  attitude  and  magnetic  torque  rods  on 
each  axis.  The  current  design  of  NPSAT-1  uses  all  of  these  sensors  to  provide  rate 
information  for  damping  and  stability  to  the  control  system  that  requires  a  complicated 
attitude  control  design.  The  purpose  of  this  attitude  control  design  simulation  is  to 
investigate  and  propose  a  control  law  utilizing  a  single  pitch  momentum  wheel  and  three 
magnetic  torque  rods.  A  further  proposal  is  to  utilize  a  constant  speed  momentum  wheel 
to  avoid  momentum  damping  and  over  speed,  replace  the  pitch  control  with  magnetic 
torquers,  and  develop  a  Kalman  filter  estimator  to  provide  all  the  required  angular  rates. 

Specific  to  problems  targeted  by  the  simulations,  external  disturbance  moments 
will  cause  errors  in  the  spacecraft's  attitude.  These  errors  will  be  kept  within  the  required 
pointing  limits  if  the  attitude  control  system  is  properly  designed.  The  five  major 
disturbance  moments  worth  consideration  are  1)  Solar  Pressure,  2)  Gravity  Gradient,  3) 
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Magnetic  Moment,  4)  Aerodynamic,  and  5)  Internal  Inertia  and  Torque.  For  the  purposes 
of  this  paper,  both  magnetic  and  aerodynamic  moments  will  be  discounted.  For  design 
simplicity,  it  will  be  assumed  that  the  solar  pressure  moment  can  be  modeled  as  a 
constant  torque  about  each  body  axis. 

General  design  of  the  NPSAT-1  will  be  a  three-axis  stabilized,  biased  momentum 
system  implementing  four  sensors  and  two  actuators.  Because  the  spacecraft  will  be 
Earth  oriented  in  a  polar  orbit,  attitude  control  about  all  three  axes  is  determined  to  be 
within  one  degree  (+/-1.00).  The  attitude  knowledge,  however,  must  be  within  0.1 
degrees  (+/-0.10)  with  jitter  minimized  to  less  than  one  milli-radian  at  1-kHz.  Navigation 
will  be  accomplished  onboard  via  a  GPS  receiver,  and  external  spacecraft  guidance  is  not 
anticipated. 

NPSAT-1  is  designed  with  four  information-gathering  sensors  and  two  actuators 
to  guide  and  control  its  movement.  These  sensors  are:  3-axis  magnetometer,  3-axis  star 
sensor,  earth  sensors,  and  a  GPS  receiver.  These  four  sensors  are  defined  to  gather 
information  to  feed  the  information  stream  supporting  the  Kalman  filter  and  the  ADCS. 

In  addition  to  the  four  sensors,  NPSAT-1  is  designed  with  two  actuators  to  guide 
and  control  its  movement.  These  are  the  Pitch  Momentum  Wheel  and  the  Magnetic 
Torquers.  A  single  momentum  wheel  will  be  placed  along  the  pitch  axis  of  the  NPSAT- 
1 .  Angular  momentum  generated  by  the  spinning  wheel  will  provide  gyroscopic  stability 
to  the  spacecraft  about  the  roll  and  yaw  axes.  Pitch  will  be  controlled  by  spinning  the 
wheel  up  or  down.  Roll  and  yaw  as  well  as  momentum  dumping  is  to  be  controlled  by 
the  torque  rods  via  control  laws.  A  set  of  three  torquers,  one  on  each  axis,  will  be 
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employed  to  initially  orient  the  NPSAT-1.  They  also  are  to  be  used  to  spin  up  the  pitch 
momentum  wheel  for  momentum  dumping  and  to  control  roll  and  yaw  via  control  laws. 

The  results  obtained  in  this  thesis  are  quite  extraordinary.  The  controller  uses  a 
magnetic  torque  actuator  to  create  the  required  torques.  The  linear  principle  of 
superposition  allowed  removal  of  wheel  speed  changing,  creating  a  constant  speed  wheel. 
The  system  was  well  behaved.  The  orbit  inclination  is  also  a  concern.  This  approach 
will  probably  have  problems  with  equatorial  or  polar  orbits. 

This  thesis  shows  that  a  properly  designed  optimal  rate  estimator  Kalman  filter  is 
effective  and  able  to  estimate  body  rate  from  a  single  star  sensor.  In  addition,  the  initial 
results  prove  that  a  single  sensor  coupled  with  a  proper  rate  estimator  design  can  be  used 
as  a  backup,  or  even  primary,  attitude  determination  process. 

For  NPSAT-1,  a  single  star  sensor  estimator  will  be  an  addition  to  the  control 
system.  This  approach  is  necessary,  especially  during  the  initial  launch,  where  after 
initial  launch  the  spacecraft  will  be  tumbling  at  some  pitch,  roll,  and  yaw  rates.  The 
magnetometer  and  the  magnetic  torquers  control,  but  would  require,  an  additional  vector 
which  the  star  sensor  can  provide. 
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I. 


INTRODUCTION 


The  micro-satellite  space  industry  has  grown  at  an  alanning  rate  over  the  last 
decade  and  will  continue  to  develop  into  an  advanced  and  sophisticated  space 
marketplace.  The  majority  of  small  satellites  in  orbit  or  under  development  were 
designed  and  manufactured  by  universities  in  collaboration  with  Aerospace  Corporation. 

To  maximize  revenue,  today's  satellites  are  designed  to  have  extended  lifetimes 
and  precise  attitude  control  systems.  The  majority  of  satellites  in  orbit  is  exceeding  their 
life  cycles  and  operating  beyond  their  designed  lifetime.  However,  some  of  the  satellite 
components  are  developing  problems  directly  related  to  the  degradation  of  their  critical 
components,  such  as  the  gyroscopic  rate  sensors.  These  rate  sensors  are  a  critical  element 
that  is  vital  to  overall  functioning  and  provides  an  attitude  determination  and  control 
system  (ADCS)  to  the  spacecraft.  To  counter  this  problem,  rate  estimator  development  is 
needed  from  different  sensor  sources  that  fuse  multiple  data  streams  together  for  use  as 
an  attitude  determination  and  control  system. 

Many  factors  can  shorten  the  lifetime  of  a  satellite,  such  as:  space  environment, 
disturbance  torques,  orbit  type,  and  spacecraft  complexity.  A  spacecraft  becomes 
complex  when  deployment  mechanisms  and  moving  parts  are  incorporated  into  its 
design.  A  rate  gyroscope,  for  example,  is  a  constantly  spinning  piece  of  equipment, 
which  provides  rate  infonnation  to  the  attitude  detennination  and  control  system.  By 
way  of  example,  if  a  satellite  is  designed  to  have  an  operational  life  of  15  years,  the 
reliability  of  its  rate  gyroscope  decreases  over  the  operational  lifetime,  more  so  than  a 
non-moving  or  solid-state  piece  of  hardware.  Rate  infonnation  provides  robustness, 
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damping,  and  stability  to  a  control  system.  If  this  stability  is  lost,  a  secular  disturbance 
torque  could  conceivably  destabilize  the  satellite  [1]. 

In  the  project  that  this  paper  is  based  on,  students  at  the  Naval  Postgraduate 
School  designed  a  micro-satellite  functioning  as  a  space-borne  platform  for  multi-spectral 
imaging.  NPS  Space  Systems  Engineering  students  designed  the  satellite  for  the  Naval 
Postgraduate  Satellite  program  (NPSAT-1)  as  a  research,  design,  fabrication,  testing, 
integration,  launch,  and  on-orbit  implementation  of  a  Low  Earth  Orbiting  (LEO)  satellite. 
The  satellite’s  mission  was  conceptualized  in  support  of  a  science  payload  dedicated  to 
multi-spectral  imaging  of  the  earth's  Aurora. 

This  paper  will  discuss  in  some  detail  the  design  of  two  of  the  main  components 
of  the  NPSAT-1,  namely,  the  Attitude  Dynamics  Control  System  and  Kalman  Filter  Rate 
Estimator.  The  discussion  of  these  components  will  include  their  conceptualization, 
theoretical  origins,  tests  and  simulations,  a  presentation  of  the  results,  conclusions  drawn 
from  the  results  and  how  they  apply  to  the  final  design  of  NPSAT-1,  and,  suggestions  for 
future  research. 

A.  SATELLITE  FAILURE  ISSUES  AND  IN-TIME  SOLUTIONS 

Many  satellites  have  been  developed  since  Explorer  I  was  launched  in  1959, 
which  perform  a  multitude  of  missions.  Mission  possibilities  include:  communications, 
mapping,  and  weather  observation.  Some  of  these  satellites  are  spin-stabilized,  dual  spin- 
stabilized,  and  three-axis  stabilized.  Due  to  the  growing  need  for  power,  most  of  today's 
spacecraft  are  three-axis  stabilized  and  will  be  the  type  of  spacecraft  considered  in  this 
study.  Typically,  communication  satellites  are  in  a  near  geo-stationary  orbit  for  middle- 
to-low  latitude  coverage,  and  in  a  Molniya  orbit  for  high  latitude  coverage. 
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As  we  move  into  the  21st  century,  we  become  increasingly  dependent  on  satellites. 
If  a  spacecraft  fails  before  its  design  life  runs  out,  the  resulting  loss  in  revenue  can  be 
staggering.  Although  many  factors  can  contribute  to  spacecraft  failure,  the  most  common 
limiting  factors  are  fuel,  batteries,  and  solar  cell  degradation.  This,  of  course,  is  highly 
dependent  upon  the  altitude  and  mission  of  the  satellite.  For  long  duration  missions,  it  is 
not  uncommon  for  a  rate  gyroscope  to  fail.  This  can  be  a  critical  failure  if  the  satellite's 
attitude  control  laws  require  inputs  from  this  device. 

Operational  information  and  data  received  from  a  satellite  are  generally  ignored 
until  a  critical  outage  occurs.  As  an  example,  the  unexpected  and  premature  failure  of  the 
Galaxy  IV  spacecraft  temporarily  left  millions  of  people  without  pager  service.  This 
failure  translates  into  a  multi-million  dollar  loss  of  financial  revenue  for  businesses 
highly  dependent  on  this  technology.  That  is  why  when  a  satellite  fails  and  its  life 
expectancy  is  threatened,  every  effort  will  be  made  to  save  it. 

Recently,  the  SOHO  spacecraft  lost  its  orientation  after  suffering  from  multiple 
gyroscope  failures.  With  that,  engineers  were  forced  to  create  a  software  package  that 
would  override  the  failed  hardware.  It  took  nearly  six  months  to  write  and  test  the  code 
before  the  control  system  of  the  billion-dollar  satellite  was  restored  as  a  successful 
uplink.  SOHO  is  now  able  to  autonomously  maintain  proper  attitude  relative  to  the  sun 
using  its  star  tracker  as  the  primary  control  sensor.  SOHO  is  testimony  that  after  a 
certain  hardware  failure  there  will  be  a  great  need  to  develop  software  upgrades  to  current 
systems  on  orbit  in  order  to  extend  its  life  cycle. 

A  part  of  this  research  is  to  devise  a  reliable  method  of  obtaining  angular  rates 
from  a  star  tracker  in  the  event  of  a  gyroscope  failure.  Since  star  sensors  only  measure 
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errors  in  angle,  rates  will  have  to  be  derived  based  on  these  measurements.  At  first 
glance,  it  might  seem  reasonable  to  simply  take  the  derivative  of  the  angles  to  get  the 
required  rate  infonnation,  but  doing  so  would  only  amplify  the  noise  effects  of  the  sensor. 
Pseudo-rate  modulators  can  be  used  to  derive  rates,  but  the  accuracy  of  these  modulators 
is  strictly  a  function  of  sensor  noise.  A  Kalman  filter  however,  can  determine  angular 
rates  as  well  as  reduce  noise,  which  is  inherent  in  any  sensor. 

Spacecraft  attitude  detennination  algorithms  traditionally  relied  only  on  rate 
gyros.  These  gyros  are  highly  reliable  but  deteriorate  over  time  and  degrade  system 
performance  dramatically  as  evident  in  the  complete  replacement  of  the  gyros  onboard 
the  Hubble  Space  Telescope.  Recent  advances  in  star  trackers  have  allowed  several 
research  projects  to  develop  Kalman  filter  algorithm  simulations.  Such  research  utilizes 
the  single  most  accurate  instrument  on  spacecraft  thus  far,  the  star  sensor. 

B.  OVERVIEW  OF  NPSAT-1  ATTITUDE  DYNAMICS  CONTROL  SYSTEM 

DESIGN 

The  design  of  the  Attitude  Control  System  for  NPSAT-1  was  developed  by 
students  during  the  preliminary  design  phase  of  AA4871  Spacecraft  Design  II 
Engineering,  with  the  assistance  of  Professor  Barry  Leonard.  During  this  phase,  certain 
design  criteria  were  specified  based  on  the  mission  requirements  of  the  spacecraft.  The 
proposed  control  law,  modeling,  and  simulation  using  MATLAB  were  developed  during 
AA3818  Spacecraft  Attitude,  Dynamics,  and  Controls,  as  a  design  project  with  the 
assistance  of  Professor  Barry  Leonard.  The  proposed  MATLAB's  SIMULINK  block 
diagram  of  the  attitude  control  system  is  shown  below  in  Ligure  1. 
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Figure  1:  Attitude  Control  System  (Systems  Design  Diagram) 

The  goal  of  the  project  described  in  this  paper  required  the  proposal  of  a  control 
law  that  can  substitute  magnetic  torquers  that  use  the  Earth’s  magnetic  field  in  place  of 
reaction  jets. 

The  following  steps  were  required: 

•  Add  the  magnetic  field  in  orbital  coordinates  to  the  simulation. 

•  Translate  the  field  components  into  body  coordinates. 

•  Assume  “ideal”  magnetic  torquers  mx,  mv,  and  mz  followed  by  a  saturation 
characteristic. 

•  Form  the  cross  product  of  m  =  (mx,  my,  m2f  with  the  B  body  coordinates 
to  generate  the  control  torques  Mc  =  (Mxc  MyCi  M:c)r. 

•  Disconnect  the  reaction  jet  control  torques  from  the  spacecraft  and  replace 
it  with  the  magnetic  control  torques. 

•  Derive  the  magnetic  torquers  with  the  following  control  laws: 

o  mx  =  0 

o  my  =  0 

o  m:  =  0 
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•  Remove  the  disturbance  torques  M  =  (Mdx,  Mdy,  Mdz)r  and  demonstrate 
acquisition  from  6  =  5°,  (j)  =  5°,  y/  =  5°.  If  possible,  acquisitions  from 
large  angles  and  perturbation  torques  should  be  perfonned. 

C.  DESIGN  OF  KALMAN  FILTER  RATE  ESTIMATOR 

In  1960,  R.  E.  Kalman  published  his  famous  paper  describing  a  recursive  solution 
to  the  discrete-data  linear  filtering  problem.  Since  that  time,  due  in  large  part  to  advances 
in  digital  computing,  the  Kalman  filter  has  been  the  subject  of  extensive  research  and 
application,  particularly  in  the  area  of  autonomous  or  assisted  navigation. 

In  this  paper,  the  Kalman  filter  will  be  applied  to  a  rate  estimator  that  is  integral  to 
controlling  a  small  satellite.  The  Kalman  filter  is  an  important  component  in  satellite 
navigation  because  it  has  the  potential  to  determine  angular  rates  as  well  as  reduce  noise, 
which  is  inherent  in  any  sensor.  A  major  part  of  this  project,  and  as  demonstrated  in  this 
paper,  is  the  simulation  of  the  Kalman  filter  and  its  interaction  with  the  navigational 
components  of  the  small  satellite  [2]. 

D.  LITERATURE  REFERENCE  NOTE 

The  primary  literature  source  for  this  paper  is  the  “NPSAT-1  Preliminary  Design 
Report”  dated  September  15,  1999  [1]  and  noted  in  the  list  of  references  at  the  end  of  this 
paper.  Because  of  the  pervasive  use  of  this  document  throughout  the  project,  the  author 
does  not  explicitly  cite  each  reference  to  this  report. 
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II.  SPACE  ENVIRONMENT 


The  space  environment  consists  of  all  things  related  to  being  in  space,  such  as 
radiation,  the  effects  of  the  sun,  and  degradation  of  the  satellite  system  itself.  A 
combination  of  the  satellite’s  mission  and  the  space  environment  itself  aid  in  determining 
the  satellite’s  orbit  type.  Key  factors  in  selecting  a  satellite’s  orbit  type  is  its  altitude 
selection  based  on  the  satellite’s  radiation  environment  and  mission  purpose.  More 
precisely,  the  satellite’s  radiation  environment  undergoes  a  substantial  change  at  about 
1000  kilometers.  Below  this  altitude,  the  atmosphere  will  quickly  clear  out  charges 
particles  so  the  radiation  density  is  relatively  low.  Above  this  altitude  are  the  Van  Allen 
belts  of  trapped  radiation  that  can  significantly  reduce  the  lifetime  of  components. 

Mission  orbits  can  therefore  be  separated  into  two  types:  low  earth  orbits  (LEO) 
below  1000  km.  and  geosynchronous  orbits  (GEO)  that  are  above  1000  km.  Because  the 
NPSAT-1  is  designed  to  study  the  Earth’s  aurora,  a  low  earth  orbit  appears  to  be  the  most 
practical  choice. 

Orbit  selection  is  a  trading  process  balancing  the  satellite’s  mission  assignment 
against  such  factors  as  payload  size  and  weight,  launch  cost,  and  design  lifetime.  Once 
we  have  selected  the  LEO  orbit  category  as  being  the  most  reasonable,  the  orbit  itself  can 
be  refined  further.  There  are  two  possibilities  in  orbit  style  that  may  be  attractive:  a 
standard  LEO  that  is  almost  circular,  and  a  Molniya  orbit  that  is  semi-synchronous  and 
eccentric.  The  advantage  of  an  eccentric  orbit  is  that  at  apogee  the  satellite’s  velocity  is 
lower  and  offers  more  time  at  apogee.  Since  we  are  studying  the  Earth’s  aurora,  this 
might  be  an  attractive  solution,  as  the  satellite  would  have  more  time  to  record  data  if 
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positioned  at  apogee  above  the  geographical  area  of  the  aurora.  These  are  the  two  orbit 
solutions  covered  in  this  research.  For  purposes  of  clarity,  both  orbits  are  earth 
referenced. 

A.  DISTURBANCE  TORQUES 

External  disturbance  moments  will  cause  errors  in  the  spacecraft's  attitude.  These 
errors  will  be  kept  within  the  required  pointing  limits  if  the  attitude  control  system  is 
properly  designed.  The  five  major  disturbance  moments  worth  consideration  are  1)  Solar 
Pressure,  2)  Gravity  Gradient,  3)  Magnetic  Moment,  4)  Aerodynamic,  and  5)  Internal 
Inertia  and  Torque.  For  the  purposes  of  this  paper,  both  magnetic  and  aerodynamic 
moments  will  be  discounted.  For  design  simplicity,  it  will  be  assumed  that  the  solar 
pressure  moment  can  be  modeled  as  a  constant  torque  about  each  body  axis.  These 
moments  were  found  to  be  [2], 

Tggx=  3QV(a 

Tgg},=m2d(rz-ix),  (1) 

^=o- 

Note  that  the  symbols  in  Equation  (1)  as  well  as  all  other  symbols  used  in  this 
paper  are  defined  in  the  section  LIST  OF  SYMBOLS  detailed  in  the  front-matter  of  this 
paper. 

1,  Modeling  Disturbance  Torques 

The  following  disturbance  torques  was  estimated  using  SMAD.  The  total  worst- 
case  disturbance  torque  was  used  to  size  the  actuator  components.  Additionally,  a 
MATLAB  and  SIMULINK  model  of  the  ADCS  sub-system  that  was  generated  took  into 
account  all  of  the  disturbances  discussed  in  the  following  subsections. 
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a)  Magnetic 

The  Magnetic  Disturbance  Torque  is  cyclic  throughout  orbit.  The 
magnetic  torque  is  altitude  dependent  with  higher  torque  at  lower  altitudes.  The  torque 
can  be  estimated  using  the  following  equation, 

T.  =  DB  =  d(2M/r,),  (2) 


where, 


2 

D  =  Residual  Dipole  of  the  Spacecraft  =  2  A-m~ 

B  =  Earth’s  Magnetic  Field  ~  2M/R3  (for  a  near-polar  orbit) 

M=  Earth’s  Magnetic  Moment  =  9.00  x  1015  Tesla-m3 
R  =  Radius  of  Orbit  =  6878  x  103  m  (Altitude  =  500  km) 

The  worst-case  Magnetic  Field  Torque  was  calculated  to  be  1.1 1  x  10'4  N-m. 

b)  Aerodynamic 

The  Aerodynamic  Torque  can  be  estimated  using  the  following  equation. 
This  torque  is  altitude  dependent  and  is  at  its  worst  for  lower  altitudes, 

Ta  =  0.5 pCdAV\Cpa  -  Cm),  (3) 


where, 


12  3 

p  =  Atmospheric  Density  =  2.8  x  10'  kg/m 
Cd=  Coefficient  of  Drag  =  2.5 
A  =  Surface  Area  of  Spacecraft  =  1.0  m2 
V=  Spacecraft  Velocity  =  7613  m/s 

(■ Cpa  -  Cm)  =  Center  of  Aerodynamic  Pressure/Mass  offset  =  0.2  m 
The  worst-case  Aerodynamic  Torque  was  calculated  to  be  1.52  x  10'5  N-m. 

c)  Solar  Radiation 

Solar  Radiation  Torque  is  a  cyclic  torque  that  varies  throughout  the  orbit. 
This  torque,  however,  was  not  dependent  on  altitude.  It  was  estimated  using  the 
following  equation, 
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where, 

Fs  =  Solar  Constant  =  1399  W/m2 
c  =  Speed  of  Light  =  3  x  108  m/s 

2 

As  =  Spacecraft  Surface  Area  =1.0  m 

q  =  Reflectance  Factor  =  0.7 

i  =  Angle  of  Incidence  of  Sun  =  0  (worst  case) 

(Cps  -  Cm)  =  Center  of  Solar  Pressure/Mass  offset  =  0.2  m 

Using  the  worst-case  values  above,  we  obtain  the  Solar  Radiation  Torque  to  be  a 

maximum  of  1.59  x  10'6  N-m. 

d)  Gravity  Gradient 

In  a  Molniya  orbit,  gravity  gradient  moments  will  be  greatest  at  perigee. 
The  Gravity  Gradient  Torque  is  given  by, 

Tgg  =  | rxagdm  .  (5) 

The  gravitational  acceleration  is, 

ag  =  -GM9 

R  is  the  distance  to  the  center  of  mass  of  the  satellite  measured  from  the  center  of  the 
Earth  and  it  is  given  by, 

R0=R63.  (7) 

Equation  (7)  expressed  in  body  coordinates  is, 

Rb=C°R0  •  (8) 


For  now,  Rh  will  be  written  as, 

Rh  =  xbt  +  Yb2  +  Zb3 .  (9) 

Taking  the  cross  product,  we  obtain 
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(10) 


rxR  =  (yZ  -  zY)bx  +  (zX  -  xZ)b 2  +  ( xY  -  yX)b3 


From  the  binomial  theorem,  the  following  expression  is  obtained 


-  _  -s  1  xX  +  yY  +  zZ 

R  +  r  =  —  -3- 


(11) 


It  can  be  shown  that  the  orbital  angular  velocity  is  just 


a  = 


GM„ 


(12) 


Equations  (9),  (10),  and  (11)  can  be  substituted  into  Equation  (12)  to  get  the  following 
expression 

Tggx  =  -3Q 2my  91  xy  -  ryz] , 

Tm,  =  -3Q2  [ 9(1  x  -E)  +  <!>Lxy  +  Ixz  ] ,  (13) 

Tgg:  =3Q  \<t>lxz+91yz). 


These  three  equations  were  derived  using  small  angle  approximations  and  ignoring 
second  order  terms. 

B.  ATTITUDE  DYNAMICS  CONTROL  SYSTEM  (ADCS)  GUIDANCE  AND 

CONTROL  CONCEPT 

The  NPS  Aurora  Satellite  (NPSAT-1)  requires  a  scheme  to  sense  disturbances  and 
a  method  of  reacting  to  disturbance.  This  section  discusses  the  general  design  concept 
for  the  ADCS. 

1,  General 

The  NPSAT-1  will  be  a  three-axis  stabilized,  biased  momentum  system 
implementing  four  sensors  and  two  actuators.  Because  the  spacecraft  will  be  Earth 
oriented  in  a  polar  orbit,  attitude  control  about  all  three  axes  is  detennined  to  be  within 
one  degree  (+/-1.00).  The  attitude  knowledge,  however,  must  be  within  0.1  degrees  (+/- 
0.1°),  with  jitter  minimized  to  less  than  one  milli-radian  at  1-kHz.  Navigation  will  be 
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accomplished  onboard  via  a  GPS  receiver,  and  external  spacecraft  guidance  is  not 
anticipated. 

2.  Total  Worst-Case  Disturbance 

The  worst-case  external  disturbance  torque  is  calculated  as  the  sum  of  the  above 
external  disturbance  torques.  While  this  method  may  result  in  an  over  design  of  the 
ADCS,  it  accounts  for  the  low  probability  event  of  ah  worst-case  external  torques 
occurring  simultaneously  and  in  the  same  direction.  Therefore,  the  total  worst-case 
disturbance  torque  for  NPSAT-1  is, 

Td  =  Tm  +  Ta  +  Tsp  =  1.04  x  1(T4  N-m.  (14) 

C.  COMPONENTS 

NPSAT-1  is  designed  with  four,  infonnation-gathering  sensors  and  two  actuators 
to  guide  and  control  its  movement.  This  section  discusses  each  of  the  components 
required  by  the  design  concept  for  information  gathering  and  control. 

1.  Sensors 

There  are  four  sensors  defined  to  gather  information  for  the  NPSAT-1.  These  are 
the  Magnetometer,  Star  Sensor,  Earth  Sensors,  and  a  GPS  Receiver. 
a)  3-Axis  Magnetometer 

A  three-axis  magnetometer  will  be  used  to  detect  Earth’s  magnetic  field  to 
determine  the  NPSAT-l’s  attitude  to  approximately  5-10  degrees  accuracy.  Output 
magnetometer  sensor  will  be  used  as  the  primary  method  of  attitude  determination  after 
separation  from  the  launch  vehicle  and  initially  orient  the  NPSAT-1.  Once  initial 
orientation  is  completed,  the  Earth  sensors  will  take  over  for  final  attitude  acquisition. 
The  magnetometer  will  also  contain  inputs  to  the  magnetic  torquers,  to  be  discussed  later. 


12 


b)  3-Axis  Star  Sensor 

A  star  camera  will  be  used  to  automatically  determine  the  NPSAT-l’s 
attitude  in  three  axes.  This  will  be  used  as  a  backup  method  of  initial  acquisition.  The 
star  sensor  is  limited  to  operate  only  when  pointed  at  space  without  interference  from  the 
Sun,  Earth,  Moon,  or  other  objects.  Accuracy  of  the  star  sensor  is  well  above  the 
requirement  of  0.1  degrees  and  provides  yaw  information  not  available  from  the  Earth 
sensors.  The  star  sensor  could  provide  attitude  knowledge  a  majority  of  the  time,  in  and 
out  of  eclipse,  while  a  sun  sensor  would  not  operate. 

c)  Earth  Sensors 

Two  sets  of  static  infrared  horizon  sensors  will  provide  the  primary 
control  input  to  the  NPSAT-l’s  actuators.  The  first  set  will  be  ‘coarse’  sensors  with  a 
field  of  view  of  20°.  This  set  will  be  used  the  for  Earth  horizon  acquisition  following 
initial  orientation  guidance  from  the  magnetometers.  The  second  set  of  horizon  sensors 
are  ‘fine’  sensors  and  will  be  used  for  the  primary  control  of  the  NPSAT-l’s  attitude. 
Both  sets  of  horizon  sensors  provide  roll  and  pitch  knowledge.  The  coarse  sensors  are  set 
to  within  1.15°,  while  the  fine  sensors  are  set  to  within  0.15°  accuracy. 

d)  GPS  Receiver 

A  GPS  receiver  is  to  be  used  for  obtaining  NPSAT-1  location  information. 
The  receiver  would  aid  the  user  in  locating  the  satellite,  and  hence,  the  observables  from 
the  payloads.  The  altitude  of  the  NPSAT-1  is  low  enough  that  GPS  satellites  will  be 
accessible  to  the  GPS  receiver.  Signals  from  this  system  will  give  ground  controllers 
exact  position  and  velocity  of  the  satellite,  thus  providing  accurate  orbital  information.  If 
GPS  technology  allows,  it  might  be  possible  to  gain  backup  attitude  information  from 
differential  GPS  through  judicious  placement  of  4  to  5  GPS  antennas.  For  example,  an 
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antenna  could  be  placed  on  the  zenith  face  and  on  each  tip  of  the  unfolding  solar  array. 
Theoretically,  this  would  provide  the  necessary  spacing  for  spacecraft  attitude  to  be 
calculated. 

2.  Actuators 

In  addition  to  the  four  sensors,  NPSAT-1  is  designed  with  two  actuators  to  guide 
and  control  its  movement.  These  are  the  Pitch  Momentum  Wheel  and  the  Magnetic 
Torquers.  A  discussion  of  these  actuators  is  the  topic  of  this  section. 

a)  Pitch  Momentum  Wheel 

A  single  momentum  wheel  will  be  placed  along  the  pitch  axis  of  the 
NPSAT-1.  Angular  momentum  generated  by  the  spinning  wheel  will  provide  gyroscopic 
stability  to  the  spacecraft  about  the  roll  and  yaw  axes.  Pitch  will  be  controlled  by 
spinning  the  wheel  up  or  down.  Roll  and  yaw  as  well  as  momentum  dumping  is  to  be 
controlled  by  the  torque  rods  via  control  laws. 

b)  Magnetic  Torquers  (Torque  Rods) 

A  set  of  three  torquers,  one  on  each  axis,  will  be  employed  to  initially 
orient  the  NPSAT-1.  They  also  are  to  be  used  to  spin  up  the  pitch  momentum  wheel,  for 
momentum  dumping,  and  to  control  roll  and  yaw  via  control  laws.  The  torque  rods  are  to 
be  double  wound  for  redundancy. 

D.  ATTITUDE  CONTROL  SYSTEM  DESIGN  SPECIFICATIONS 

The  development  of  Spacecraft  Attitude,  Dynamics,  and  Controls  design  is  to 
build  a  suitable  attitude  control  law.  The  attitude  control  system  will  be  designed 
according  to  the  following  parameters. 

1.  Satellite  Specifications 

o  LEO  orbit 
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o  Two  earth  sensors  aligned  with  the  body  axes 
o  One  star  sensor 

o  Three-axis  magnetometer 

o  Global  Positioning  System  (GPS) 

o  Roll  inertia,  Ixx= 24.67  kg-m2 

o  Pitch  inertia,  Iyv= 22.63  kg-m- 

o  Yaw  inertia,  Izz=  1 1  kg-m" 

o  550  km  altitude,  circular  orbit 

o  Kalman  filter  rate  estimator  (future  development) 

2.  Other  Specifications 

o  Nadir  pointing  to  within  +/-  1° 

o  4  arc-second  noise  level  for  each  sensor 

o  Three-axis  stabilized 

o  1 -year  design  life 

3.  Assumptions 

o  Small  and  large  angle  (up  to  10  r/s)  approximations 

o  Orbital  angular  velocity  and  acceleration  known  for  each  sensor 
measurement 

o  Constant  solar  pressure  moments 

o  No  slewing  requirement 

o  Satellite  is  modeled  as  a  rigid  body 

4.  Attitude  Control  System  Design  Considerations 

In  order  to  achieve  1.0°  pointing  accuracy,  a  constant  speed  momentum  wheel, 
and  three  magnetic  torquers  whose  momentum  vectors  coincide  with  the  body  axes  will 
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be  employed.  As  disturbance  moments  cause  errors  in  attitude,  off-axis  components  of 
reaction  wheel  angular  momentum  will  cause  internal  torques  that  must  be  accounted  for. 
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III.  MODELING  AND  SIMULATION 


Many  types  of  kinematical  transfonnation  methods  are  in  use  in  various  types  of 
research,  but  the  most  popular  are:  direction  cosine  matrices  (DCM),  Euler  angles,  and 
quaternion.  Quaternion  is  popular  since  it  involves  only  a  single  rotation  about  an  Eigen- 
axis.  On  the  other  hand,  making  small  angle  approximations  and  setting  second  order 
tenns  to  zero,  DCM's  are  easily  employed  and  are  used  in  this  analysis.  Transformations 
from  one  frame  to  another  are  perfonned  to  facilitate  calculations.  For  example,  the 
latitude  and  longitude  of  stars  in  the  star  catalog  have  all  been  programmed  within  a 
celestial  frame,  but  measurements  will  be  made  in  the  body  frame.  Therefore,  proper 
attitude  detennination  relies  on  a  simple  transfonnation  [3]. 

In  the  field  of  attitude  control,  it  is  often  required  to  express  an  inertial  quantity  as 
a  body  frame  quantity.  For  example,  the  inertial  angular  velocity  derived  from  the  Euler 
moment  equations  must  be  expressed  in  body  coordinates,  and  then  integrated  to  get  the 
Euler  angles.  The  body  frame,  orbital  frame,  and  inertial  frame  are  the  three  reference 
frames  are  used  in  the  derivation  of  equations  of  motion:.  The  origin  of  these  three 
frames  will  all  be  located  at  the  spacecraft's  center  of  mass.  In  the  right-hand  set,  the 
orbital  reference  frame,  the  Z-axis  points  at  the  center  of  the  Earth,  the  X-axis  points  in 
the  satellite's  direction  of  motion,  and  the  Y-axis  is  normal  to  the  orbital  plane, 
completing.  In  the  left-hand  set,  the  body  reference  frame  is  attached  to  the  spacecraft; 
therefore,  the  Euler  angles  represent  the  deviation  of  the  body  reference  frame  from  the 
orbital  reference  frame.  On-board  sensors  measure  these  Euler  angles.  The  inertial 
frame  remains  fixed  in  Earth  space  such  that  the  inertial  Y-axis  coincides  with  the  orbital 
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Y-axis.  The  celestial  frame  is  an  additional  reference  frame  alluded  to  earlier.  The  Z- 


axis  of  this  frame  points  north  and  the  X-axis  points  in  the  direction  of  the  vernal 
equinox.  Although  the  X-axis  precesses  (albeit  very  slowly),  the  assumption  is  that  it  is 
fixed  in  space. 

A.  STAR  SENSOR  CHARACTERISTICS 

The  star  sensor  model  used  in  this  simulation  was  designed  in  accordance  with  the 
specifications  outlined  in  Chapter  II.  Table  1  summarizes  the  characteristics  of  the  star 
sensor  and,  for  completeness;  additional  assumptions  have  been  made  that  will  be 
consistent  with  current  technology. 


Technology 

Charged  Coupled  Device  (CCD) 

FOV 

10°  x 10° 

Accuracy 

~10  arcsec 

#  Stars  in  Catalog 

4000 

Sampling  Rate 

0.1  Hz  (current  technology  is  faster) 

Noise 

4  arcsec  (magnitude=6) 

Solar  Exclusion  Angle 

30°  w/sun  shade 

Table  1:  Star  Tracker  Characteristics. 


The  noise  level  shown  in  Table  1  is  inherent  to  the  star  tracker  itself  and  it  is  treated  as  a 
zero-mean  Gaussian  white  sequence. 

B.  GPS  SYSTEMS 

Three-axis  attitude  determination  requires  two  separate  line-of-sights  (LOS)  with 
angular  separation  near  90°  for  increased  accuracy.  In  this  simulation,  the  optical  axis  of 
each  star  sensor  will  be  aligned  with  the  body  axes,  and  at  each  discrete  time  step  a  star 
sensor  will  be  selected  at  random  for  attitude  determination.  Although  this  sequence  of 
events  pennits  only  one  LOS  per  time  step,  attitude  is  readily  detennined  over 
consecutive  time  steps.  Since  only  one  star  will  be  in  the  sensor  FOV  at  any  particular 
time,  measurements  can  only  be  made  about  two  axes. 
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c. 


REFERENCE  FRAMES 


The  simulation  discussed  in  this  section  makes  use  of  three  previously  discussed 
frames  of  reference;  the  body  frame,  the  inertial  frame,  and  the  orbital  (LVLH)  frame.  A 
fourth  frame  of  reference,  the  earth  frame,  is  introduced  to  complete  the  simulation. 
Orientation  of  the  frames  is  as  follows;  the  body  frame  is  fixed  to  the  spacecraft  and 
aligned  with  the  satellite's  principal  moments  of  inertia,  the  inertial  frame  is  celestially 
fixed,  while  the  angular  velocities  of  the  orbit  and  earth  frames  with  respect  to  the  inertial 


The  angular  velocity  of  the  body  frame  with  respect  to  the  orbital  frame  is  given  by 

o  (b  b=  +  [6C"  ]9/f  +  [AC°  ]/7o3  ,  (16) 

where  n  is  just  an  intermediate  frame.  Since  on-board  sensors  make  measurements  with 

respect  to  the  body-frame,  all  of  the  above  rotation  rates  will  be  transfonned  to  the  body- 

frame.  As  can  be  seen  from  Equations  (15)  and  (16),  the  C  matrices  (DCM),  perform  this 

transfonnation.  For  example,  the  DCM  that  transforms  the  orbital  frame  to  the  body 

frame  is  given  by  the  following  3-2-1  rotations 

c9cy/  cdsif v  -s9 

bC°=  -  c<fsy/  +  stjisOcy/  c  </>c  y/  +  s  0s  9s  y/  s0c9  .  (12) 

s0sy/  +  c0s9cy/  -  s 0c y/  +  c (j>s 9s y/  c9c(f> 

Expanding  Equation  (16)  and  getting  rid  of  2nd  order  terms, 

0  co  b=  +  yss<fc9^jt>2  +  ^c9c y/  -  .  (18) 

From  Equation  (18),  the  Euler  rates  can  now  be  determined, 
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(19) 


(f)  =  +  ®2  tan#sin(zi  +  ®3  tan^cos^  , 

8  =  ®2  cos^  -  ®3  sin  $  , 

sintfi  +  costfi 

w  =  — - 1 - - - —  . 

cos  8 

These  rates  can  be  integrated  to  give  the  Euler  angles.  These  angles  represent  the 
deviation  of  the  body  frame  with  respect  to  the  orbit  frame.  The  fine  horizon  sensor  will 
detect  these  errors  and  feed  them  to  the  momentum  wheel  that,  in  turn,  will  act  as  an 
actuator  to  correct  this  error.  The  magnetometer  and  the  star  sensor  will  act  in  a  similar 
fashion;  however,  the  star  sensor  will  not  be  part  of  the  control  loop.  Instead,  its  attitude 
information  will  be  sent  down  to  the  Alaska  ground  station  for  attitude  knowledge.  The 
operational  purpose  of  this  attitude  control  system  is  to  keep  the  body  frame  aligned  with 
the  orbital  frame.  The  following  table  lists  the  characteristics  of  this  3-axis  stabilized 
micro-satellite. 


Altitude 

550  km 

Period 

5677.0  s 

na 

0.0011068  rad/s 

^xx 

22.222  kg-m2 

4v 

21.387  kg-m2 

4, 

17.056  kg-m2 

Pointing  Knowledge 

.1°  all  axes 

Pointing  Accuracy 

1°  all  axes 

Table  2:  Satellite  Properties 

D.  SATELLITE  DYNAMICS 

Since  it  is  assumed  that  the  spacecraft  rotates  about  its  principal  moments  of 
inertia,  the  spacecraft's  angular  momentum  is  given  by  [6] 


H  = 


/ 

0 

0  0 


0  0 


]yy  0 


CO 


(20) 
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Notice  that  to  use  Equation  (20),  the  angular  velocity  must  be  inertial  and  it  must  be 
expressed  in  body  coordinates;  using  small  angle  approximations,  it  can  be  determined  by 
adding  Equation  (15)  and  (18) 

lmb  ={^>-Q.0Y]pi+{d-D.0%  +{¥  +  ^J%-  (21) 

It  is  also  important  to  note  that  H,  in  Equation  (20),  is  the  total  angular  momentum  of  the 
satellite  in  expressed  in  body  coordinates,  including  the  momentum  wheel.  It  can  be 
broken  up  as  follows 

H  =  Hb  +  hw .  (22) 

Each  angular  momentum  component  rotates  about  its  own  center  of  mass,  and  the 
momentum  wheel  only  rotates  about  the  negative  pitch  axis,  so  that  hw  =  -hb2 .  The 
inertial  rate  of  change  of  angular  momentum  is  just  equal  to  the  external  moment 

M  =  —H  =  —H+‘abxH  .  (23) 

dt  dt 

The  external  moment,  M,  in  Equation  (23)  represents  the  sum  of  all  of  the  external 
moments,  including  gravity  gradient,  aerodynamic,  solar  radiation  pressure,  magnetic 
control,  and  pitch  wheel  control.  The  gravity  gradient  moment  is  given  as  Equation  (16) 

Mgg  =  3  n2ju=  -  />,  +  3  n20e{ia  -  ixx)b2 .  (24) 

The  other  moments  will  be  derived  later.  Substituting  Equations  (21),  (22),  and  (24)  into 
Equation  (23),  the  following  results  are  obtained 

Mdx  =  1tJ>  +  [~  4^o {hz  -  Iyy)+  ^■oh\l>  +  [“  ^o([xx  ~  1  yy  +  Iz:)+  ’ 

Mdy  =  Ij - 3 n2„ (/zz  - Ixx p  +  h,  (25) 

Mdz  =  hzV  +  [nl(lyy  ~  Ixx)+  ohV  +  [Q0(7«  “  1  yy  +  1  zz)~  ^  • 

As  can  be  seen  from  Equation  (20),  the  pitch  moment  equation  is  independent  of  roll  and 
yaw.  These  moment  equations  are  in  agreement  with  Equation  ( 1 5). 
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E. 


PITCH  CONTROL 


The  solution  to  Equation  (25)  is  oscillatory,  if  h  =0,  i.e.  no  wheel  control,  the 
following  occurs 


Figure  2:  Pitch  response  with  no  damping. 


The  response  in  Figure  2  can  be  dampened  out  by  controlling  the  momentum  wheel,  h  , 
according  to  the  following  control  law 


h  =  kgTgd  +  kgO  .  (26) 

Substituting  Equation  (26)  into  Equation  (25)  and  taking  the  Laplace  transform,  the 

following  transfer  function  is  derived 


®fr) 

MdyiS) 


s  + 


kdTe 


yy 


s  +  - 


(27) 


The  plot  in  Figure  2  assumes  an  initial  condition  of  9=  5°,  but  Equation  (27)  is  based  on 
the  initial  conditions  being  zero.  It  is  of  no  consequence,  however,  since  the  initial 
conditions  will  not  affect  the  characteristic  equation.  The  steady  state  error  can  be 
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determined  for  a  step  input  disturbance  torque  using  the  final  value  theorem 
Equation  (19)  and  it  can  be  shown  that 


= 


M 


dy 


2  ‘ 


IyyWe 


(28) 


The  allowable  steady  state  error  given  by  the  requirements  in  Table  2  is  ±1°,  but  .8°  will 
be  chosen  as  the  design  margin.  In  addition,  the  constants  in  Equation  (41)  will  be 
determined  based  on  a  critically  damped  system.  The  following  table  lists  the  pitch 
control  properties. 


steady  state  error,  8SS 

±.8° 

worse  case  pitch  disturbance  moment,  Mdv 

-.0001  N-m 

damping  ratio, ^ 

1 

natural  frequency,  roe 

0.0182995  rad/s 

pitch  time  constant,  te 

109.0031  s 

pitch  auto-pilot  gain,  k0 

0.00718096  N-m/rad 

Table  3:  Pitch  Controller  Properties 

Using  the  controller  gains  given  in  Table  3,  the  following  response  is  obtained 


Pitch  Response  to  -1  e-4  step  input,  with  controller 

4 

3 

1  2 

1 

0 

-1. 

\ 

\ 

5  100  200  300  400  500  600  700  800  900  1000 

t 

Figure  3:  Pitch  response  with  damping. 
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F.  ROLL- YAW  RESPONSE 

The  horizon  sensor  is  capable  of  detecting  both  pitch  and  roll  errors,  so  both  of 
these  angles  will  be  available  for  state  feedback.  The  momentum  wheel  gives  the 
spacecraft  gyroscopic  stiffness  along  the  pitch  axis,  so  that  when  a  roll  or  a  yaw  error 
occurs,  the  result  is  mutation  about  the  pitch  axis,  similar  to  a  spinning  top.  Assuming 
that  /z»max(/Q0),  it  can  be  shown  that  Equations  (25)  can  be  simplified  to 


Mdx  =  Jxx<i  +  hnot  +  hy/, 

Mdz  =  IZZW  +  hQ.0y/  -hij) . 


(29) 


Taking  the  Laplace  transform  of  Equations  (29)  the  result  is  four  transfer  functions; 
however,  the  design  transfer  functions  are  the  following 


ZM  .  4a2+^V> 

Md:{s)  [lxxs2  +  noh\l:zs2  +  noh)+  h2ss  ’ 


$(■?)  _ 4a!  +  aoh 

Md-As )  (/VA2  +  noh\lz:s2  +  Q0/?)+  h2ss 

The  yaw  and  roll  response  of  these  two  transfer  functions  is  shown  in  Figure  4. 


(30) 
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Figure  4:  Undamped  roll  and  yaw  responses  to  step  input. 

The  roll  response  is  within  the  pointing  requirements,  but  the  yaw  response 
exceeds  1°  after  40  minutes.  This  is  assuming  that  the  initial  conditions  for  roll  and  yaw 
are  both  zero.  If  the  angular  momentum  of  the  wheel  is  increased,  both  responses  would 
be  within  limits;  for  example,  if  h=15Nms,  operationally,  no  roll-yaw  damping  would  be 
required.  For  acquisition,  however,  roll  and  yaw  damping  will  be  required  since  there 
will  be  an  initial  error  upon  launch  vehicle  separation.  It  is  interesting  to  note  that  the 
period  of  the  responses  is  simply  the  period  of  the  orbit;  there  is  also  a  short  period 
response,  which  is  related  to  momentum  wheel  precession.  Evidence  of  roll-yaw 
coupling  can  also  be  seen  since  the  two  plots  are  90°  out  of  phase. 


G.  ATTITUDE  CONTROL  DETERMINATION 

Many  types  of  control  laws  are  available,  which  can  conceivably  satisfy  this 
satellite's  pointing  requirements.  Some  common  control  laws  are:  1)  proportional,  2) 
proportional  plus  derivative,  3)  proportional  plus  integral  plus  derivative,  and  4)  optimal. 
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Each  of  these  controllers  has  its  own  unique  characteristics;  however,  as  long  as  the 
controllers  maintain  proper  spacecraft  attitude,  controllers  that  are  more  exotic  will  not  be 
required.  In  fact,  it  will  be  shown  that  the  gains  of  a  simple  PD  controller  can  be  adjusted 
to  minimize  overshoot  and  settling  time, 

hx  =kvJ  +  kx</>, 

hy=kvye  +  kye,  (31) 

K  =  kvzijr  +  kzy/  . 


These  control  laws  are  expressed  as  the  rate  of  change  of  reaction  wheel  angular 
momentum,  or  reaction  wheel  torque,  and  they  are  part  of  the  feedback  loop.  As  can  be 
seen,  these  internal  torque  equations  are  a  function  of  the  measured  Euler  angles  and 
rates.  If  the  resulting  set  of  equations  is  completely  de-coupled,  and  the  Laplace 
transfonn  is  taken,  the  following  result  is  obtained, 


1 

$(4 _ A _ 

4(4  2  kvx  4Q 2{i  -ix)-ahy+k/ 

s  +— —s  + - - - - - 

4  4 

l 

0(4  _ 4 _ 

4(4  2  k  3Q2(4  -Iz)  +  k  ’ 

S  + — —  S  + - — 

4  4 

j_ 

T(4  _ 4 _ 

4(4  2  kvz  n2(-ix+iy)-nhv+kz  ' 

5  + — —  S  + - ; ; - 

4  4 


(32) 


For  this  particular  analysis,  it  is  assumed  that  the  orbital  angular  velocity  is  locally 
constant.  The  objective  is  to  determine  suitable  position  and  rate  feedback  gains  that  will 
increase  spacecraft  robustness.  The  nominal  characteristic  equation  for  any  second  order 
system  has  the  following  form 

A(s)  =  s2  +2a>ngs  +  a>l .  (33) 
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The  natural  frequency  is  denoted  as  con  and  ^  is  the  damping  factor,  which  will 
be  chosen  to  be  one.  Each  of  the  denominators  in  Equation  (32)  will  be  equated  to 
Equation  (33).  Solving  for  the  coefficients,  the  result  is  two  equations  and  three 
unknowns.  The  third  equation  makes  use  of  the  final  value  theorem,  and  it  is  given  by 
the  following 

/ (oo)  =  lim  / ( t )  =  lim  sF(s) .  (34) 

t—>  oo  s  — >0 


The  pointing  requirements  for  this  satellite  require  a  steady  state  pointing  accuracy  of 
0.1°  about  each  axis.  By  applying  the  final  value  theorem  to  Equation  (32)  and  if  the 
external  disturbance  torques  can  be  approximated  as  a  step  input,  position  feedback  gains 
can  be  detennined  from  the  following  equations 


Tx  -4Q 2(Iy  -Ix)</>ss  +OhJss 

K 

T  -3Q2(/x  -Iz)0ss 
ky=- - a - ’ 

°ss 

Tz  -n2(-fx  + 1 y  )¥ss  +  D.hyy/ss 


(35) 


The  'ss'  subscript  denotes  steady  state  and  the  design  torques  represent  a  worst- 
case  scenario.  It  can  be  seen  from  Equation  (35)  that  the  position  feedback  gains  are  not 
constant;  they  will  vary  as  a  function  of  orbital  position.  The  natural  frequency  for  roll, 
pitch,  and  yaw  can  now  be  determined  by  taking  the  square  root  of  the  last  term  in  the 
denominator.  Once  this  is  found,  the  velocity  feedback  gains  can  be  calculated  from  the 
following  expressions 

kvx  —  2a>nxIx , 

kyy  =  2 (0„yly  ,  (36) 

kyz  —  z  • 
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In  a  similar  manner  to  the  position  feedback  gains,  the  velocity  feedback  gains  also  vary 
with  time.  Equation  (36)  and  Figure  7  each  depicts  the  time  varying  nature  of  the  PD 
gains  over  one  orbit,  specifically  during  perigee.  As  expected,  the  pitch  and  pitch  rate 
gains  are  much  higher  than  the  roll  and  yaw  gains  [7], 

H.  KINEMATICS 

The  inertial  frame,  orbital  frame,  and  body  frame  are  the  three  reference  frames 
are  used  in  the  derivation  of  equations  of  motion.  Direction  cosine  matrices  (DCM)  are 
used  to  transform  between  coordinate  systems.  These  matrices  are  given  by  Equations 
(15)  and  (16)  above.  The  orbital  frame  of  reference  is  oriented  such  that  the  x-axis  points 
in  the  direction  of  the  velocity  vector,  the  z-axis  points  towards  the  center  of  the  Earth, 
and  the  y-axis  completes  the  right-hand  set.  The  body  frame  is  aligned  with  the  orbital 
frame,  as  that  is  the  direction  of  motion.  The  following  3-2-1  transformations  from  the 
orbital  frame  to  the  body  frame  are  given  by  Equation  (17).  The  orbital  frame  rotates  at  a 
rate  of  Q(?)  with  respect  to  the  inertial  reference  frame,  or 


' cb°  =  -Qo2  .  (37) 

To  perform  angular  momentum  calculations,  the  inertial  angular  velocity  is  expressed  in 
body  coordinates.  This  is  represented  by 


i  -\b  i  -\o  .  o  -\b 
CO  =  CO  +  CO  . 


(38) 


Expressed  in  body  coordinates,  the  angular  velocity  of  the  orbital  frame  with  respect  to 
the  inertial  frame,  is 


ld)°h  =-bc°nd2 .  (39) 

Angular  velocity  of  the  body  frame  with  respect  to  the  orbital  frame,  expressed  in  body 
coordinates,  is 
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°abb  =  #1  +  C(^)C(e)0h2+bCoVo 3 . 


(40) 


The  /72unit  vector  belongs  to  an  intermediate  reference  frame.  If  Equations  (39)  and  (40) 
are  substituted  into  (38),  the  following  result  is  obtained 

lwh  =  {<j> -Q.y/)bx  +  {6-D.)b2  +  {¥  +  Q0)b3 .  (41) 

Equation  (41)  is  a  simplified  expression  using  small  angle  approximations  and  neglecting 
second  order  terms  [6],  [7],  [8], 

I.  DERIVATION  OF  EQUATIONS  OF  MOTION 

In  detennining  the  attitude  of  a  satellite,  the  common  approach  is  to  translate  all  factors 
into  the  body  coordinate  system,  since  on-board  sensors  are  designed  to  detect  errors  with 
respect  to  the  body  frame.  As  noted  previously,  Equation  (34)  was  obtained  using  small 
angle  approximations  where  errors  are  represented  by:  (f>  =  roll  error,  0  =  pitch  error,  and 

i//  =  yaw  error.  Total  spacecraft  angular  momentum  can  be  separated  into  two  angular 
momentum  vectors  for  the  spacecraft  body  and  the  reaction  wheels,  which  is  given  by  the 
expression 

H  =  Hb+Hw.  (42) 

If  it  is  assumed  that  cross  products  of  inertia  are  negligible,  then 

Hb=l‘a)h.  (43) 

Note  that  when  calculating  the  angular  momentum  of  the  satellite  about  its  center  of 
mass,  inertial  angular  rates  must  be  used  rather  than  body  rates.  Substituting  Equations 
(41)  and  (42)  into  Equation  (43),  the  total  spacecraft  angular  momentum  is 

H  =  (ij-  Ixn  ¥ + hx  %  +  (iye- iyn  +  hy  %  +  {iz¥+ i.a <t> + hz  )b3 .  (44) 

Next,  the  relation 
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(45) 


d  '  -  d  h  -  ,  _h  - 
dt  dt 

is  used  to  determine  the  Euler  moment  equations.  Neglecting  second  order  terms  and 
gravity  gradient  moments,  Euler  equations  for  this  spacecraft  are 

r  =  ij + 4Q2  (/,  - /.  V  -  n hj  -  n h: + n(- /,  +  /  .  -  />  -  //.»/> + ho  -  ip  ys + K , 

Ty  =  1,6  +  3Q2  {lx  -L)9  +  hxxj/  +  Qhi//  +  ClhJ  -  hj  -  Ip  +  hy ,  (46) 

Tz  =  Lys  +  tE  (■ -I  X+Iy  V  -  Qhyy/  +  Qhx  +  fl(lx  -  Iy  +  !,}/>-  hx6  +  hj  +  Ip(j>  +  h, . 

These  equations  describe  the  motion  of  the  spacecraft  when  subject  to  external 
disturbance  torques.  Rates  of  change  of  angular  momentum  of  each  reaction  wheel  will 
be  used  to  counteract  the  disturbance  moments,  maintaining  the  required  pointing 
accuracy.  However,  solving  for  the  Euler  angles  is  non-trivial  as  all  three  differential 
equations  are  second  order  and  coupled.  If  the  cross  products  of  inertia  are  not 
negligible,  the  equations  of  motion  become  [6],  [7],  [8], 

Ta  =Tx+{G-a-3Q.26)Ixy+{0.2¥  +  yI  +  a^Ixz+(c-2aG-2D.jlyz, 

Tyi  =Ty  +(-2Q^  +  2  n2<f,  +  0-ny,)Ixy  +3Cl2Ixz  +(2<p-Cl2y/  +  yI  +  Q.</>)Iyz,  (47) 

T  ,  =T  +{2Q.G-Q.2)I  +(-2Q2^  +  (zS-q^)/  +(6> -Q-3Q2 G)I  . 

z  i  z  xy  xz  yz 

J.  CONTROL  LAWS 

To  determine  control  laws  for  the  magnetic  torquers,  one  must  define  what  needs 
to  be  controlled.  For  the  satellite,  the  best  approach  would  be  to  join  all  the  equations  in 
a  matrix  fonn  and  only  then  linearize  them.  An  engineer  would  come  up  with  a  MIMO 
control  law  that  would  control  the  system  as  a  whole  and  not  by  parts.  This  approach  has 
one  important  advantage:  it  does  not  matter  if  the  state  variables  are  coupled. 

The  only  concern  becomes  the  linearization.  As  a  side  effect,  one  can  use  the 
model  of  the  plant  in  a  Kalman  filter  configured  as  an  observer  or,  if  the  noise  sources  are 
not  relevant,  use  a  deterministic  observer.  The  engineer  could  then  use  the  wealth  of 
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methods  that  are  present  in  the  linear  state  space  approach,  especially  adaptive,  robust, 
and  optimal,  control  techniques. 

The  control  law  proposed  here  is  by  no  means  the  best  one,  only  suitable  for  the 
model.  Control  laws  were  developed  for  pitch,  by  using  the  rate  of  momentum  of  a 
wheel,  and  for  roll/yaw,  by  using  reaction  jets.  In  turn,  both  the  reaction  jets  and  the 
momentum  wheel  laws,  generated  torque  commands.  These  were  fed  into  the  physical 
devices  that  would  produce  the  actual  torque.  It  is  reasonable  to  argue  that  any  physical 
device  that  can  produce  the  toque  required  by  the  control  law  would  be  adequate.  Let  us 
keep  this  reasoning  in  mind  when  going  over  the  next  paragraphs. 

Suppose  that  it  was  possible  to  generate,  with  the  magnetic  torquers,  the  same 
torque  that  is  generated  by  the  reaction  jets.  Imagine  that  there  is  a  box  that  can  accept 
the  torque  required  and  outputs  the  commands  for  the  torquers.  One  could  then  connect 
the  output  of  the  reaction  jets  control  law  (i.e.,  a  number  that  requests  the  amount  of 
torque  required)  and  feed  it  into  that  magic  box.  Now,  we  can  bring  this  forward  a  step  if 
the  control  law  assumes  the  system  is  linear.  Of  several  implications,  one  of  the  most 
important  is  part  of  the  very  definition  of  linearity:  the  superposition  principle  must 
apply.  Next,  we  look  at  the  pitch  control.  The  pitch  controller  generates  a  number  that 
represents  the  amount  of  torque  the  pitch  wheel  must  generate.  If  we  take  this  number 
and  sum  it  up  with  the  output  by  the  reaction  jets  controller,  we  would  have  the  total 
torque  needed  to  control  the  spacecraft  attitude. 

If  we  take  the  total  torque  that  is  needed  to  control  the  spacecraft  and  feed  it  to  the 
controller,  the  torquers  could  produce  all  the  torque  needed.  No  change  in  the  angular 
velocity  of  the  wheel  would  be  needed,  and  no  more  momentum  dumping,  and  no  more 
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wheel  over  speed  and  a  larger  MTTF.  However,  open  questions  remain,  such  as;  what 
goes  inside  the  control  algorithm,  and  are  the  torquers  able  to  generate  all  the  torque 
needed  to  control  the  spacecraft? 

K.  MEAN  TIME  TO  FAILURE 

1.  Using  the  Earth’s  Magnetic  Field 

Suppose  that  one  has  the  Earth  magnetic  field  vector  B  and  the  amount  of  torque 
desired  T </.  We  need  to  know  the  value  m  needed  for  the  magnetic  torquers  that  will 
produce  Tp  in  the  presence  of  B.  Note  that  the  torque  produced  may  be,  as  will  be  seen, 
different  from  the  desired  one 


Tp  =  m  x  B  (48) 

From  Equation  (48)  follows  Equation  (49).  This  equation  will  give  the  direction 
of  the  magnetic  torquers  vector.  However,  the  engineer  must  be  cautious:  T d  will  only  be 
the  output  of  the  torquers  when  T,/  _L  B.  In  all  the  other  circumstances,  the  best  torque 
obtainable  is  the  portion  of  T*  which  is  perpendicular  to  B.  That  portion  of  T j  is  then 
called  T p 


m  =  BxT(. 


(49) 


The  amplitude  of  m  can  be  found  using  Equation  (50) 


m 


B 


sin  6  ||b|| 

mB 


(50) 


The  engineer  will  then  have  to  make  a  choice:  is  T d  replaced  in  Equation  (50) 
with  Tp?  The  question  can  only  be  answered  with  further  analysis  and  simulations.  In 
this  work,  Equation  (50)  was  used  with  a  slight  variation;  shown  in  Equation  (51),  note 
that  the  Tp  was  changed  to  T c/, 
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(51) 


This  law  leaves  a  question  open:  what  happens  when  the  Earth’s  magnetic  field  is 
parallel  to  Tf?  Nothing,  because  no  torque  can  be  generated,  resulting  in  lack  of  control 
over  the  spacecraft.  However,  as  the  time  progress,  the  magnetic  field  moves  with 
respect  to  the  spacecraft  and  control  is  regained.  How  the  shortage  is  expected  to  be 
remains  an  open  question.  What  are  the  consequences  of  T d  ^  T P1  The  answer  to  this 
question  depends  on  more  analysis  and  simulations,  ft  seems  that  for  somewhat  low 
precision  (0.1°)  the  difference  probably  can  be  compensated  by  the  controller,  fn  the 
simulations  the  steady  state  for  small  and  large  angles  was  e  =  (0  0  0  l)r.  The  results 
with  disturbances  were  similar  to  the  one  produced  with  the  reaction  jets.  B/,  does  have 
discontinuities.  What  happens  with  the  control?  Every  time  the  B  measured  in  the  body 
coordinates  changes  its  signal,  the  control  law  would  need  to  have  discontinuities  also. 
Since  the  discontinuities  are  when  the  T d  is  almost  parallel  to  B,  the  best  solution  is  to 
turn  off  the  controller  and  let  the  spacecraft  drift  until  this  condition  vanishes.  This  is 
probably  the  best  approach  because  of  the  magnetic  field  sensors  noise  and  miss 
alignment. 

There  is  another  concern  in  the  actual  implementation:  the  torquers  cannot  be 
turned  on  all  the  time.  The  torquers  must  be  turned  off  in  order  to  perform  accurate 
magnetic  field  measurements.  If  the  torquers  are  active  only  half  of  the  time,  then  the 
torquers  commands  must  be  multiplied  by  a  factor  of  two.  ft  is  easy  to  come  out  with  the 
correction  factors  for  other  operation  conditions. 
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L.  DISCRETE  KALMAN  FILTER 

The  Kalman  Filter  is  a  recursive  algorithm  for  estimating  a  state  vector  given  past 
estimates  and  current  measurements  with  noise.  With  a  system  model  of  the  plant 
dynamics  and  sensor  noise,  the  filter  will  minimize  the  mean  square  error.  Since  it  was 
first  development  in  1960  by  R.  E.  Kalman,  the  filter  has  been  used  in  numerous  fields  of 
study  and  many  sources  exist  that  walk  through  the  derivation  of  his  work.  For 
simplicity,  only  the  resulting  equations  will  be  shown  here  [3],  [8]. 

The  filter  itself  is  a  two  step  process,  a  prediction  followed  by  an  update.  In  this 
simulation,  a  single  measurement  is  used  as  the  initial  prediction.  The  simulation  code 
therefore  first  computes  the  Kalman  Gain  with  the  initial  covariance  prediction  before 
updating  the  covariance  prediction  and  then  making  a  new  prediction. 

1.  The  Process  to  Be  Estimated 

The  Kalman  filter  addresses  the  general  problem  of  trying  to  estimate  the  state 
x  e  9? N  of  a  discrete -time  controlled  process  that  is  governed  by  the  linear  stochastic 
difference  equation 

xk+i  =  Axk  +  Blh  +  wk »  (52) 

with  a  measurement  zk  e  9? N  that  is 

zk=Hkxk+vk-  (53) 

The  random  variables  wk  and  vk  represent  the  process  and  measurement  noise 
respectively.  They  are  assumed  independent  of  each  other,  white,  and  with  nonnal 
probability  distributions 
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p(w)~  N{Q,Q), 
p[v)- Nib,!!). 


(54) 


The  n  x  n  matrix  A  in  the  difference  Equation  (68)  relates  the  state  at  time  step  k 
to  the  state  at  step  k+1,  in  the  absence  of  either  a  driving  function  or  process  noise.  The 
matrix  B  relates  the  control  input  to  the  state  x.  The  matrix  H  in  the  measurement 
Equation  (53)  relates  the  state  to  the  measurement  z*. 

2.  The  Computational  Origins  of  the  Filter 

We  define  xk  e  iH N  (note  the  “super  minus”)  to  be  our  a  priori  state  estimate  at 

step  k  given  knowledge  of  the  process  prior  to  step  k,  and  xk  e  tH5'  to  be  our  a  posteriori 

state  estimate  at  step  k  given  measurement  z/(.  We  can  then  define  a  priori  and  a 
posteriori  estimate  errors  as 

ck=xk-xk,  (55) 

and 

ck=xk-xk.  (56) 

The  a  priori  estimate  error  covariance  is  then 

Pi=E[vlT\  (57) 

and  the  a  posteriori  estimate  error  covariance  is 

Pk  =4vfi-  (58) 

In  deriving  the  equations  for  the  Kalman  filter,  we  begin  with  the  goal  of  finding 
an  equation  that  computes  an  a  posteriori  state  estimate  xk  as  a  linear  combination  of  an 
a  priori  estimate  xk  and  a  weighted  difference  between  an  actual  measurement  and  a 
measurement  prediction  as  shown  below  in  Equation  (59).  Some  justification  for 
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Equation  (59)  is  given  in  the  sub-section  “The  Probabilistic  Origins  of  the  Filter,”  found 
below 

**  =  **  +  K(zk  ~Hk*k)-  (59) 

The  difference  [zk  -  H kxk )  in  Equation  (59)  is  called  the  measurement 
innovation,  or  the  residual.  The  residual  reflects  the  discrepancy  between  the  predicted 
measurement  Hkxk  and  the  actual  measurement  zk.  A  residual  of  zero  means  that  the 
two  are  in  complete  agreement. 

The  n  x  m  matrix  K  in  Equation  (59)  is  chosen  to  be  the  gain  or  blending  factor 
that  minimizes  the  a  posteriori  error  covariance,  Equation  (58).  This  minimization  can  be 
accomplished  by  first  substituting  Equation  (59)  into  the  above  definition  for  ek , 
substituting  that  into  Equation  (58),  performing  the  indicated  expectations,  taking  the 
derivative  of  the  trace  of  the  result  with  respect  to  K,  setting  that  result  equal  to  zero,  and 
then  solving  for  K.  One  form  of  the  resulting  K  that  minimizes  Equation  (59)  is  given  by 

K„  =Pk-H;{HtPt~Hl  +  R„y 

Pi  Hi  ■  (60) 

h„p;hI  +  r„ 

Looking  at  Equation  60,  we  see  that  as  the  measurement  error  covariance  Rk 
approaches  zero,  the  gain  K  weights  the  residual  more  heavily.  Specifically 

lim  Kk  =Hk  .  (61) 

o 

On  the  other  hand,  as  the  a  priori  estimate  error  covariance  Rk  approaches  zero, 
the  gain  K  weights  the  residual  less  heavily.  Specifically 

lim  K,  =  0  .  (62) 

p;->  o 
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Another  way  of  thinking  about  the  weighting  by  K  is  that  as  the  measurement 
error  covariance  Rk  approaches  zero,  the  actual  measurement  zk  is  “trusted” 

increasingly,  while  the  predicted  measurement  Hkxk  is  trusted  increasingly  less.  On  the 
other  hand,  as  the  a  priori  error  covariance  estimate  Pk  approaches  zero,  the  actual 
measurement  zk  is  also  trusted  increasingly  less,  while  the  predicted  measurement  Hkxk 
is  increasingly  trusted. 

3.  The  Probabilistic  Origins  of  the  Filter 

The  justification  for  Equation  (59)  is  rooted  in  the  probability  of  the  a  priori 
estimate  xk  conditioned  on  all  prior  measurements  zk  (Bayes’  rule).  For  now  let  it 

suffice  to  point  out  that  the  Kalman  filter  maintains  the  first  two  moments  of  the  state 
distribution 

E[xk 

E\*k  -xk\xk  -xj 

The  a  posteriori  state  estimate  equation  60  reflects  the  mean  (the  first  moment)  of 
the  state  distribution — it  is  normally  distributed  if  the  conditions  of  Equations  (54)  are 
met.  The  a  posteriori  error  covariance  estimate,  Equation  (58),  reflects  the  variance  of 
the  state  distribution  (i.e.,  the  second  non-central  moment).  In  other  words 

p(xk  \zk)~ N[E[xk ],E\xk  - xk \xk  -  xk )r  J)  =  N(xk ,Pk).  (64) 

4.  The  Discrete  Kalman  Filter  Algorithm 

We  begin  this  section  with  a  broad  overview  covering  the  “high-level”  operation 
of  one  form  of  the  discrete  Kalman  filter.  After  presenting  this  high-level  view,  we  will 
narrow  the  focus  to  the  specific  equations  and  their  use  in  this  version  of  the  filter. 
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The  Kalman  filter  estimates  a  process  by  using  a  form  of  feedback  control:  the 
filter  estimates  the  process  state  at  some  time  and  then  obtains  feedback  in  the  fonn  of 
(noisy)  measurements.  As  such,  the  equations  for  the  Kalman  filter  fall  into  two  groups: 
1)  time  update  equations  and,  2)  measurement  update  equations.  The  time  update 
equations  are  responsible  for  projecting  forward  in  time  the  current  state  and  error 
covariance  estimates  to  obtain  the  a  priori  estimates  for  the  next  time  step.  The 
measurement  update  equations  are  responsible  for  the  feedback,  (i.e.)  for  incorporating  a 
new  measurement  into  the  a  priori  estimate  to  obtain  an  improved  a  posteriori  estimate. 

The  time  update  equations  can  also  be  thought  of  as  predictor  equations,  while  the 
measurement  update  equations  can  be  thought  of  as  corrector  equations.  Indeed  the  final 
estimation  algorithm  resembles  that  of  a  predictor-corrector  algorithm  for  solving 
numerical  problems  as  shown  below  in  Figure  5. 


Time  Update  Measurement  Update 

(“Predict”)  (“Correct”) 


Figure  5:  The  ongoing  discrete  Kalman  filter  cycle. 

The  time  update  projects  the  current  state  estimate  ahead  in  time,  while  the  measurement 
update  adjusts  the  projected  estimate  by  an  actual  measurement  at  that  time. 

The  specific  equations  for  the  time  updates  are  presented  below  in  Equations  (65) 
while  the  measurement  updates  are  presented  below  in  Equations  (66). 
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(65) 


vk+l 


—  ^k^k 


+  Bu 


k  ’ 


Pk+ 1  =  AkPk  K  +Qk- 

Notice  how  the  time  update  in  Equations  (65)  project  the  state  and  covariance 
estimates  from  time  step  k  to  step  k+1 .  Ak  and  B  are  from  Equation  (52),  while  Qk  is 
from  Equation  (54).  Initial  conditions  for  the  filter  were  discussed  in  earlier  references 

=  K  +  k(z,  -  Hti;  ),  (66) 


The  first  task  during  the  measurement  update  is  to  compute  the  Kalman  gain,  Kk . 
Note  that  the  equation  given  here,  as  Equation  (66),  is  the  same  as  Equation  (60).  The 
next  step  is  to  actually  measure  the  process  to  obtain  zk,  and  then  to  generate  an  a 

posteriori  state  estimate  by  incorporating  the  measurement  as  in  Equation  (66).  Again, 
Equation  (66)  is  simply  Equation  (59)  repeated  here  for  completeness.  The  final  step  is 
to  obtain  an  a  posteriori  error  covariance  estimate  via  Equation  (66). 

After  each  time  and  measurement  update  pair,  the  process  is  repeated  with  the 
previous  a  posteriori  estimates  used  to  project  or  predict  the  new  a  priori  estimates.  This 
recursive  nature  is  one  of  the  very  appealing  features  of  the  Kalman  filter  as  it  makes 
practical  implementations  much  more  feasible  than,  say  for  example,  an  implementation 
of  a  Weiner  filter  which  is  designed  to  operate  on  all  of  the  data  directly  for  each 
estimate.  Instead,  the  Kalman  filter  recursively  conditions  the  current  estimate  on  all  of 
the  past  measurements.  Figure  6  below  offers  a  complete  picture  of  the  operation  of  the 
filter,  combining  the  high-level  diagram  of  Figure  5  with  the  equations  from  Equations 
(65)  and  Equations  (66). 
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5. 


Filter  Parameters  and  Tuning 


In  the  actual  implementation  of  the  filter,  each  of  the  error  measurement 
covariance  matrices  Rk  and  the  process  noise  Qk ,  as  given  by  Equations  (54),  might  be 
measured  prior  to  operation  of  the  filter.  In  the  case  of  the  measurement  error  covariance 
Rk ,  in  particular,  this  makes  sense-because  we  need  to  be  able  to  measure  the  process, 
while  operating  the  filter;  generally,  we  should  be  able  to  take  some  off-line  sample 
measurements  in  order  to  determine  the  variance  of  the  measurement  error. 

In  the  case  of  Qk,  oftentimes  the  choice  is  less  detenninistic.  For  example,  this 

noise  source  is  often  used  to  represent  the  uncertainty  in  the  process  model  shown  in 
Equation  (52).  Sometimes  a  very  poor  model  can  be  used  simply  by  “injecting”  enough 
uncertainty  via  the  selection  of  Qk walues.  In  this  case,  one  would  hope  that 

measurements  of  the  process  would  be  reliable. 

In  either  case,  whether  or  not  we  have  a  rational  basis  for  choosing  the 
parameters,  statistically  speaking,  superior  filter  perfonnance  can  be  obtained  by  “tuning” 
the  filter  parameters  Qk  and  R, .  The  tuning  is  usually  perfonned  off-line,  frequently 
with  the  help  of  another  distinct  Kalman  filter. 
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Initial  estimates  for  xk  and  Pk 


V 


Time  Update 
(“Predict”) 

Measurement  Update 
(“Correct”) 

(1)  Project  the  state  ahead 

(1)  Compute  the  Kalman  gain 

**+i  =  Akxk  +  Buk 

K,  =PkHl{n,P,Hl  +R,Y 

(2)  Project  error  covariance  ahead 

(2)  Update  measurement  with  zk 

Pk+i  =  +  Qk 

xk  =  xk  +  K(zk  -  Hkxk  ) 

(3)  Update  the  error  covariance 

Pt  =(l-KtH„)Pt- 

Figure  6:  A  complete  picture  of  the  operation  of  the  Kalman  filter,  combining  the  high- 
level  diagram  of  Figure  5  with  the  equations  from  Equations  (65)  and  (66). 


In  closing  we  note  that  under  conditions  where  Qk  and  R,  are  constant,  both  the 
estimation  error  covariance  Pk  and  the  Kalman  gain  Kk  will  stabilize  quickly  and  then 

remain  constant  (see  the  filter  update  equations  in  Figure  6).  If  this  is  the  case,  these 
parameters  can  be  pre-computed  by  either  running  the  filter  off-line  or,  for  example,  by 
determining  the  steady-state  Pk  value. 

It  is  frequently  the  case  however  that  the  measurement  error  does  not  remain 
constant.  For  example,  when  sighting  beacons  in  our  optoelectronic  tracker  ceiling 
panels,  the  noise  in  measurements  of  nearby  beacons  will  be  smaller  than  in  beacons  that 
are  more  distant.  In  addition,  the  process  noise  Qk  is  sometimes  changed  dynamically 

during  filter  operation  in  order  to  adjust  to  different  dynamics.  As  an  example,  in  the 
case  of  tracking  the  head  of  a  user  of  a  3D  virtual  environment  we  might  reduce  the 
magnitude  of  Qk  if  the  user  seems  to  be  moving  slowly,  and  increase  the  magnitude  if 

the  dynamics  start  changing  rapidly.  In  such  a  case  Qk  can  be  used  to  model  not  only  the 
uncertainty  in  the  model,  but  also  the  uncertainty  of  the  user’s  intentions. 
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A  6-state  discrete  Kalman  filter  has  been  chosen  to  estimate  both  position  and 
rates  from  noisy  star  sensor  data.  The  Kalman  filter  that  will  be  used  in  the  simulation  is 
represented  by 


Xk+l  =®kXk+AkUk+Wk 
zk  =  Hxk  +  vk 


(67) 


The  white  sequence  wk  for  the  plant  has  a  covariance,  Q,  while  the  sensor’s  noise 
vk  has  a  covariance,  R.  Noise  from  the  star  sensor  is  affected  by  the  magnitude  of  the 

star;  a  bright  star  is  noisier  than  a  dim  star.  The  sensor  noise  covariance  is  defined  as 
follows 

Rt=E^yt}.  (68) 

M,  DERIVATION  OF  THE  Q  MATRIX 

Solving  for  the  covariance  of  the  plant  noise  is  no  trivial  matter.  In  this 
simulation,  the  Q  matrix  will  vary  with  each  time  step.  The  fonnal  definition  of  the  plant 
noise  covariance  is  given  by  [3],  [10] 

Qk=E[wkwTk\  (69) 

It  can  be  shown  that  Equation  (69)  must  satisfy  the  following  matrix  differential  equation 

Qk  =  AugQk  +  QXug  +bwbt.  (70) 

The  augmented  Aaug  matrix  is  defined  from  x  =  (A  -  BF)x  +  Bud  as  the  quantity,  A-BF, 

and  the  power  spectral  density  matrix  associated  with  the  forcing  function  u  is  denoted 
by  W. 

The  solution  to  Equation  (70)  is  greatly  simplified  for  the  time  invariant  case.  It 
proceeds  as  follows, 
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(71) 


bwbt 

bt 


By  taking  the  matrix  exponential  of  Equation  (71),  the  following  result  is  obtained 

fi=  •"  .  (72) 

_°  x 

The  upper  left  partition  can  be  neglected  for  this  analysis.  The  plant  noise  covariance 
matrix  can  now  be  determined  by  multiplying  the  upper  right  partition  of  Equation  (72) 
by  x  •  This  method  was  first  fonnulated  by  Van  Loan  in  1978. 

N.  KALMAN  ALGORITHM 

Before  entering  the  Kalman  filter  loop,  an  initial  estimate  x0 ,  and  its  error 
covariance  P0  ,  must  chosen.  The  superscript  will  represent  the  predicted  estimate 
while  the  ,A’  notation  denotes  estimation.  The  discrete  Kalman  filter  is,  in  essence,  just  a 
computer  algorithm  that  derives  optimal  estimates  from  discrete  measurements. 
Although  there  are  different  forms  of  the  discrete  Kalman  filter,  the  most  fundamental 
form  starts  with  the  Kalman  gain  calculation,  which  is  given  by  [3],  [10] 

Gk  =  PkHTk  (j HkPkHTk  +  Rk r1  .  (73) 


The  value  of  this  gain  matrix  will  vary  with  each  time  step.  Next,  it  is  required  to  update 
the  estimate  using  star  sensor  data  to  detennine  the  accuracy  of  this  new  estimate.  These 


equations  are  given  as 


xk  =xk  +  Gk  \zk  -Hkxk  j, 
Pk={l-GkHk)Pk. 


(74) 
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Equation  (75),  shown  below,  illustrates  the  recursive  nature  of  the  discrete  Kalman  filter. 
A  favorable  characteristic  of  any  recursive  filter  is  that  there  is  no  need  to  store  past 
measurements. 

Equation  (74)  is  the  actual  output  of  the  discrete  Kalman  filter.  It  estimates  both 
attitude  angles  and  attitude  rates  given  only  star  sensor  angle  information.  Not  only  does 
it  derive  rates,  but  it  also  mitigates  sensor  noise  effects.  Lastly,  it  is  necessary  to  project 
ahead  and  estimate  the  state  for  the  next  time  step.  The  predictive  equations  are  as 
follows 


*a+i  =  ®a*a  +  aa«a> 
Pk+ 1  =  T1  A  A  ®  A  +  Qk  ■ 


(75) 


It  is  interesting  to  note  that  in  Equation  (75),  the  deterministic  forcing  function 
has  been  included.  This  forcing  function  consists  of  known  reaction  wheel  moments, 
which  can  be  measured  by  the  reaction  wheel  motor  assembly.  If  this  detenninistic  term 
is  not  included,  the  rate  estimator  is  unable  to  accurately  estimate  satellite  rates  near 
perigee. 

For  the  purpose  of  analysis  and  proper  tuning,  it  is  helpful  to  look  at  the  time- 
varying  nature  of  both  the  Q  and  P  matrices  over  a  period  of  one  orbit.  Since  the  off- 
diagonal  elements  of  these  matrices  are  small,  only  the  diagonal  elements  will  be 
examined.  These  elements  are  shown  in  Figure  16. 

O.  STATE  SPACE  EQUATIONS  OF  MOTION 

The  equations  of  motion  completely  describe  the  motion  of  the  satellite  when 
subjected  to  both  internal  and  external  disturbance  moments.  If  the  body  accelerations 
are  solved  for  in  Equation  (47),  the  following  result  is  obtained  [8],  [10] 
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For  computational  reasons,  it  is  desirable  to  reduce  these  second  order  equations  to  first 
order  equations  by  making  the  following  state  variable  substitutions 

x  =  [</>  (/>  0  0  y/  y/] 7 .  (77) 


With  these  definitions,  we  can  translate  the  satellite  dynamic  equations  into  the  following 
matrix  form 


x  ~  Ax  +  Bu  . 


A  is  the  plant  matrix  and  it  is  given  by 
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B  is  the  control  matrix  given  by 

"0  0  0 

1  0  0 

0  0  0 
B  = 

0  1  0 
0  0  0 
0  0  1 


(78) 


(79) 


(80) 


u  is  going  to  be  the  input  reaction  wheel  control,  and  will  have  the  following  form 

u  =  -Fx  +  ud .  (81) 
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F  will  be  the  PD  controller  gain  matrix  and  ud  will  represent  the  summation  of  the  solar 
pressure  moments  and  the  internal  reaction  wheel  moments.  F  and  ud  are  given  by 


Substituting  Equation  (82)  into  Equation  (78),  the  following  equation  of  motion  is 
obtained 

x  =  (A-BF)x  +  Bud .  (83) 

Equation  (83)  is  equivalent  to  the  equations  of  motion. 

P.  THE  MODELING  APPROACH 

The  derivations  excerpted  in  this  section  were  employed  as  part  of  the  foundation 
for  the  simulations  in  this  project  and  were  originally  compiled  in  a  paper  by  Henry  D. 
Travis  at  the  Naval  Postgraduate  School.  For  the  sake  of  completeness,  the  full  reference 
for  this  compilation  is:  Travis,  Henry  D.,  “Attitude  Determination  Using  Star  Tracker 
Data  with  Kalman  Filters,”  Thesis,  Naval  Postgraduate  School,  December  2001. 

1.  Astrodynamics 

a)  Equations  of  Motion 

In  this  case,  the  satellite  is  traveling  along  a  Molniya  orbit  in  the  plane  described 
by  the  radius  and  velocity  vectors.  The  spacecraft’s  position  in  the  orbit  is  defined  by  its 
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distance  in  kilometers  from  the  center  of  the  earth,  r,  and  the  true  anomaly,  v,  which  is 
measured  from  perigee  in  radians.  Using  the  basic  equations  of  motion 


Fr  =  m(r  -rv2). 


(84) 


Fv  =  m(r  v  -  2 r  v)  , 

the  forces  in  both  the  radial  and  tangential  direction  can  be  determined.  Since  the  mass 
of  the  spacecraft  is  considered  a  point  mass,  the  force  per  unit  mass  can  be  written  as 

Fr  I m  =  ~nl r 2.  (85) 

Because  of  the  nature  of  a  two-body  problem  the  only  force  acting  on  the  point 
mass  is  in  the  radial  direction,  the  angular  force  F  is  zero.  This  leaves 

r  =  r v1  -  fi / r2 , 


v  =  -2  r  v  /  r  . 


(86) 


b)  Modeling  Molniya  Orbit 

We  first  model  the  system  in  the  fonn 

x  =  Ax  +  Bu  , 
y  =  Cx  +  Du  , 

by  defining  the  state  vector,  x,  as 

x  =  [rfvv]T  , 


and  the  control  vector,  u,  from  Equations  (86)  as 


u  - 


rv2  -  /dir2 
-2 rv  I r 


(87) 


(88) 


(89) 


With  all  non-linearity  included  in  the  control  laws,  the  system  coefficients  can  be 
easily  written  as 
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With  the  linear  system  coefficient  matrices  defined,  the  next  step  is  to  compute 
the  state  transition  matrix,  ®,  and  the  convolution  matrix,  A,  using  the  ‘c2d’  function  in 
MATLAB.  Thus,  the  entire  Molniya  orbit  can  be  described  using  the  discrete  equation 

At+i  =  ®xk  +  Am*  .  (91) 

This  orbital  information  is  stored  for  future  reference  with  the  pitch  controller. 

2.  The  Discrete  Kalman  Filter 
a)  Definitions 

The  Kalman  Filter  is  a  recursive  algorithm  for  estimating  a  state  vector 
given  past  estimates  and  current  measurements  with  noise.  With  a  system  model  of  the 
plant  dynamics  and  sensor  noise,  the  filter  will  minimize  the  mean  square  error.  Since  it 
was  first  development  in  1960  by  R.E.  Kalman,  the  filter  has  been  used  in  numerous 
fields  of  study  and  many  sources  exist  that  walk  through  the  derivation  of  his  work.  For 
simplicity,  only  the  resulting  equations  will  be  shown  here.  Following  modem 
convention,  the  following  definitions  and  notation  will  be  used: 
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xk  k_]  =  State  vector  estimate  at  time  k  given  measurements  up  to  time  k- 1 
xk  k  =  State  vector  estimate  at  time  k  given  measurements  up  to  time  k 
e  =  x  -  x,  i,  =  State  estimate  error 

k\k 

P,  .  .  =  Prediction  of  the  covariance  of  the  state  vector 

k\k-\ 

Pk,k  =  Update  of  the  covariance  of  the  state  vector 

G  =  Kalman  gain 

Q  =  Plant  covariance  matrix 

R  =  Measurement  covariance  matrix 

O  =  State  transition  matrix 

H  =  Observation  matrix 

w  =  Zero  mean  white  Gaussian  plant  noise 

v  =  Zero  mean  white  Gaussian  measurement  noise 

z  =  Noisy  measurement 

b)  System  Model 

Before  we  can  build  the  filter,  a  system  must  be  developed  that  will 
adequately  describe  the  behavior  of  the  spacecraft.  For  simplicity,  we  have  linearized  the 
attitude  dynamics  equations  of  motion  and  used  the  small  angle  approximation.  Looking 
at  the  equations  of  angular  velocity  in  a  rotating  frame  with  a  y/  — >  9  — >  (/)  transformation 

COBi  —  G>br  "*■  ^ it m  ■>  (92) 

where 
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(94) 
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Thus,  the  time  rate  of  change  of  co  follows  as 
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(97) 


Using  the  preceding  dynamics  equations,  we  further  assume  that  the  second  time 
derivatives  of  the  angles  are  small  enough  to  be  ignored.  This  gives  the  following  linear, 
constant  coefficient  matrices  for  the  first  system,  which  deals  with  only  roll  and  yaw. 
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The  final  step  is  to  compute  the  state  transition  matrix,  ®  and  the  convolution 
integral,  A  using  the  MATLAB  command 

[  <D  ,  A  ]  =  c2d(A,B,dt);  (99) 

which  are  used  to  promulgate  the  state  vector  using  the  equation 

xk+ l=®xk+Auk.  (100) 

c)  Controllability 

This  linear  time-invariant  system  is  considered  controllable  if  an  input, 
u(t)  will  transfer  the  initial  state  of  the  system  x(0)  to  the  origin,  x(tf)=0  with  tf  finite. 
Setting  the  state  equation  from  the  previous  section  to  zero,  it  can  be  shown  that  the 
system  is  controllable  if  it  satisfies  the  condition  that  the  controllability  matrix, 

Cm=[B  AB\,  (101) 


has  an  inverse.  Computing  Cm  for  our  system,  it  is  seen  that  this  condition  is  satisfied. 
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d)  Control  Gains 

Two  methods  are  presented  for  computing  the  control  gain.  The  first  is 
the  pole  placement  method,  and  the  second  is  a  Linear  Quadratic  Regulator  method,  or 
LQR. 

Pole  placement  method  has  been  relied  on  for  decades  as  a  means  of  ensuring  the 
poles  of  the  closed-loop  system  are  at  desirable  locations.  Ackermann  developed  a 
procedure  for  computing  the  control  gain  and  for  single  input  systems;  this  algorithm  is 
performed  using  the  ‘acker’  command  in  MATLAB.  However,  for  a  multiple  input  case, 
the  MATLAB  ‘place’  command  is  used  instead.  The  inputs  for  the  place  command  are 
the  state  transition  matrix,  the  convolution  matrix,  and  the  desired  eigenvalues. 

Using  a  Linear-quadratic  regulator  design  for  discrete -time  systems  is  another 
method  of  computing  the  control  gain.  The  gains  from  this  method  are  considered 
optimal  since  the  state-feedback  law  u[n]  =  -kx[n]  minimizes  the  cost  function 

'=Z(  x'Qx  +  u'Ru  +  2x'  Nu),  (102) 

subject  to  the  state  dynamics 

x[n+l]  =  ®x[n]  +  Au[n],  (103) 

The  matrix  N  represents  a  relationship  between  the  system  noise,  Q,  and  the 
measurement  noise,  R,  and  is  set  to  zero  for  our  system.  Also  returned  are  the  Riccati 
equation  solution  and  the  closed-loop  eigenvalues. 

e)  Filter  Model 

We  begin  developing  the  Kalman  Filter  by  modeling  the  random  process 

**+i  =®*** +w*>  (104) 

The  process  is  observed  at  discrete  points  in  time  by  the  following  relation 
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(105) 


Since  pitch  is  decoupled  from  roll  and  yaw,  we  omit  pitch  and  pitch  rate  from  the 
first  state  vector.  The  state  vector,  x,  is  then 

x  =  ] ,  (no) 

A  second  filter  is  therefore  necessary  for  pitch  and  the  corresponding  state  vector  is 
defined  here  as 

x  p  =  [00],  (111) 
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Similarly,  the  plant  covariance  follows  as 


Q„  =  cil  x 


dt4  / 4  dt2/l 
dt2  / 2  dt 2 


(112) 


and  the  measurement  covariance  is  reduced  to 

R  =  [ste2].  (113) 

f)  Algorithm 

The  fdter  itself  is  a  two-step  process,  a  prediction  followed  by  an  update. 
In  this  simulation,  a  single  measurement  is  used  as  the  initial  prediction.  The  code 
therefore  first  computes  the  Kalman  Gain  with  the  initial  covariance  prediction  before 
updating  the  covariance  prediction  and  then  making  a  new  prediction 

G  =  Pk]k_1HT(HPk]k_lHT+Ryl, 

^k\k  —  (I  —  GH)Pk\k_]  ,  (114) 

Vi=®VDr+0- 

Similarly,  the  initial  state  vector  is  first  updated  with  the  new  Kalman  Gain  before 
a  new  prediction  of  the  state  vector  is  made  based  on  the  measurement  residual 

K\k  =xk\k-x+G(z-zk\k-x)-  (115) 

The  control,  u,  is  then  detennined  using  the  Optimal  Control  Law 

u=-k(x-xk]k_]),  (116) 


before  updating  the  state  vector  and  predicting  the  next  state  estimate  using  the  state 
transition  matrix,  O,  and  convolution  integral,  A, 

xk+l  =  Ox  +  An  , 


*k\k-l  ~  ^^k\k  +  A U  . 


(117) 
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The  same  procedure  is  followed  to  update  and  predict  the  pitch  and  pitch  rate 
estimates  with  the  new  state  vector,  xp. 

There  are  several  ways  to  initialize  the  Kalman  Filter,  including  using  an  assumed 
state,  a  measurement,  or  a  batch-processed  state.  An  assumed  state  would  be  used  either 
when  the  state  is  generally  known  without  measurement  or  when  the  error  is  permitted  to 
be  large  because  there  is  sufficient  time  for  the  larger  transient.  A  measurement  approach 
is  used  when  the  error  is  expected  to  be  small  to  begin  with.  A  batch  processed  method 
would  be  used  when  the  dynamics  of  the  system  would  cause  large  changes  from  one 
measurement  to  the  next.  The  expected  changes  will  dictate  the  number  of  measurement 
processed  to  initialize  the  filter.  Due  to  the  accuracy  of  the  star  trackers  and  assumed  low 
rate  of  change  of  the  Euler  Angles,  the  Kalman  Filter  could  use  the  first  measurement  of 
the  star  tracker  to  initialize  the  filter.  To  measure  the  responsiveness  of  our  filter,  we  will 
assume  a  zero  state  initially. 

g)  Modifications 

The  Kalman  Filter  does  an  excellent  job  of  seeing  through  the  noise  to 
provide  reliable  state  estimates.  By  taking  a  closer  look  at  the  filter,  we  see  some 
variables  that  can  be  adjusted.  Most  of  these  adjustments  will  involve  the  system 
covariance  matrix.  In  covariance  manipulation,  system  covariance  is  an  assumption  of 
how  much  noise  exists  in  the  system.  Small  covariance  correlates  to  a  small  amount  of 
noise  and  will  result  in  a  better  estimate.  If  the  system  is  subjected  to  a  large  noise  from 
an  unexpected  source,  the  filter  will  be  unable  to  track  the  transient  and  the  estimate  will 
deteriorate.  The  desired  covariance  would  be  the  smallest  possible  while  still  tracking 
any  expected  transient.  Different  variables  and  methods  of  adjusting  covariance  are: 
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reset  threshold,  discretized  residual  bias,  alpha-beta,  star  gap  response,  and  glitch/bias 
rejection. 
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IV.  RESULTS 


A.  NPSAT-1  SIMULATION  RESULTS 

Results  of  the  NPSAT-1  project  described  in  this  work,  are  presented  in  this 
section.  No  noise  analysis  was  made  due  to  the  amount  of  effort  required.  The  actual 
control  law  is  very  dependent  upon  the  cross  product  of  the  magnetic  field  vector  and  the 
amount  of  torque  desired.  Therefore,  noise  in  the  angle  between  both  vectors  is  going  to 
have  an  important  effect  on  the  results.  The  evaluation  of  this  control  law  will  not  be 
complete  without  a  study  of  the  noise  effect.  Additionally,  two  different  results  are 
presented  here:  with  magnetic  pitch  control  and  without  magnetic  pitch  control. 

1.  Small  Angles 

In  this  subsection,  the  results  of  the  system  with  small  initial  angles  and  no  noise 
are  presented.  The  system  was  initialized  with  the  quaternion  set  of 
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The  results  are  presented  in  the  following  graphs.  Figures  7,  8,  9,  and  10  display 
simulated  results  with  no  disturbance 
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Angle  Idegiee) 


E  u  lei  .sngts-s  —  No  noise,  small  angles,  i  ol  Ipric  h  yaw,  no  peiluibalion 


Figure  7:  Small  angles,  mechanic  pitch,  with  perturbation. 


Figure  8:  Small  angles,  mechanic  pitch,  with  perturbation,  expanded. 


Eulei  a  ng  les  —  No  no  ise  .  sma  II  a  nej  les,  ■  oil  p  Mch'yaw.  wH  h  peilu  ibai  io  n 


Figure  9:  Small  angles,  mechanic  pitch,  with  perturbation,  transient  response. 


58 


Figure  10:  Small  angles,  mechanic  pitch,  with  perturbation,  steady  state  response. 


Figures  11,  12,  13,  and  14  display  simulated  results  with  disturbance: 


Eulei  angles  —  1-Jo  noise,  small  angles,  i  oil  .yaw,  no  peiluibalion 


Figure  11:  Small  angles,  mechanic  pitch,  with  perturbation,  transient  response. 
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Figure  12:  Small  angles,  mechanic  pitch,  with  perturbation,  transient  response. 


Eulei  angles  —  tslo  noise,  small  s  rigid’s.  10  llyaw.  with  peiluibalcn 


Figure  13:  Small  angles,  mechanic  pitch,  with  perturbation,  transient  response. 


Figure  14:  Small  angles,  mechanic  pitch,  with  perturbation,  steady  state  response. 
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2.  Large  Angles 


In  this  subsection,  the  results  of  the  system  with  small  initial  angles  and  no  noise 
are  presented.  The  system  was  initialized  with  the  quaternion  set  of 

e. 


Since  the  most  interesting  results  for  the  large  angles  acquisition  is  when  the 
system  is  under  torque  perturbations,  only  these  results  are  presented  here.  In  addition, 
the  steady  state  is  not  relevant  for  this  analysis,  since  the  results  are  the  same  as  with  the 
small  angle  acquisition.  The  results  of  the  simulation  are  presented  as  follows: 


u.o 

0.5 

0.5 

,/l-3t0.52') 


(119) 


Figure  15:  Large  angles,  magnetic  pitch,  with  perturbation. 
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Figure  16:  Large  angles,  mechanic  pitch,  with  perturbation. 


Figure  17:  Large  angles,  mechanic  pitch,  with  perturbation. 


3.  Wheel  Bias  Failure 

An  interesting  result  would  be  if  the  satellite  could  be  controlled  without  the 
momentum  bias  wheel.  Since  the  magnetic  torquers  can  control  roll,  pitch,  and  yaw  this 
would  be  a  nice  result. 

The  attempt  of  removing  the  wheel  from  the  loop  resulted  in  instability  in  yaw. 
This  was  expected,  since  the  roll/yaw  gains  were  computed  for  the  system  with  the  wheel. 
Interesting  enough,  the  controller  could  hold  very  well  in  roll  and  pitch  (even  much  better 
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than  with  the  wheel  on).  It  is  believed  that  the  controller  coefficients  can  be  recalculated 
to  yield  a  stable  system. 

This  approach  is  interesting  in  case  of  wheel  failure.  However,  it  is  not 
recommended  for  general  control.  Since  the  removal  of  the  wheel  also  removes  stiffness 
form  the  system,  any  small  torques  can  cause  large  drift  on  the  satellite  attitude  in  small 
amount  of  time.  If  that  occurs,  the  magnetic  torquers  may  not  be  able  to  control  the 
attitude  in  a  satisfactory  manner  due  to  the  control  shortages"  created  when  the  desired 
torque  is  aligned  with  the  Earth  magnetic  field. 

B.  RATE  ESTIMATOR  SIMULATION  RESULTS 

The  results  of  the  3 -axis  star  sensor  Kalman  filter  rate  estimator  are  outlined 
below.  These  simulations  prove  that  a  single  star  sensor  can  provide  the  rate  estimator 
required  to  provide  additional  vectors  for  the  NPSAT-1  initial  attitude  determination.  For 
the  Molniya  orbit  case,  the  single  axis  star  sensor  provides  estimated  yaw  rates  based  on 
the  Kalman  filter,  which  can  be  coupled  with  other  sensors  in  future  designs  to  provide  a 
more  accurate  rate  estimate,  which  in  turn,  will  be  utilized  in  the  event  of  a  rate  gyro 
failure. 

The  results  of  the  Kalman  filter  rate  estimator  simulations  are  outlined  as  follows: 
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Figure  18:  Molniya  orbit  simulation. 
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Figures  (19),  (20),  and  (21)  represent  constant  gain  Kalman  filter  with  q=0.01 
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Figure  19:  Pitch  (q=0. 01). 
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Figure  20:  Roll  (q=0.01). 
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Figure  21:  Yaw  (q=0.01). 
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Figures  (22),  (23),  and  (24)  represent  constant  gain  Kalman  filter  with  q=10.0 
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Figure  22:  Pitch  (q=10.0). 
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Figure  23:  Roll  (q=10.0). 
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Figure  24:  Yaw  (q=10.0). 
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Figures  (25),  (26),  and  (27)  represent  constant  gain  Kalman  filter  with  q=  1000.0. 


pitch  &  pitch  est  vs  time 


Figure  25:  Pitch  (q=1000.0). 
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Figure  26:  Roll  (q=1000.0). 


Figure  27:  Yaw  (q=1000.0). 
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V. 


SUMMARY  AND  CONCLUSION 


The  process  of  the  design  and  simulation  of  a  small  satellite  is  an  essential  part  of 
technological  development  required  to  explore  and  test  new  design  techniques  and 
procedures.  When  exploring  new  technologies,  it  is  essential  to  utilize  available  and 
proven  current  technologies,  and  test  these  technologies  as  backup  systems  to 
demonstrate  design  feasibility. 

Additionally,  the  design  and  development  of  a  optimal  Kalman  filter  rate 
estimator  to  perfonn  the  essential  requirement  for  attitude  control  determination  can  be 
useful  for  future  development  of  a  more  complex  rate  estimator  necessary  to  implement 
more  advanced  and  complex  control  systems. 

A.  NPSAT-1  SUMMARY 

The  results  obtained  in  this  thesis  are  quite  extraordinary.  The  controller  uses  a 
magnetic  torque  actuator  to  create  the  required  torques. 

The  linear  principle  of  superposition  also  allowed  the  removal  of  wheel  speed 
changing,  creating  a  nice  constant  speed  wheel.  The  system  was  well  behaved.  The  orbit 
inclination  is  also  a  concern  (this  approach  will  probably  have  problems  with  equatorial 
or  polar  orbits).  There  is  an  important  detail  that  must  be  mentioned:  the  controller  dead 
zone. 

B.  KALMAN  FILTER  RATE  ESTIMATOR 

This  thesis  has  shown  that  a  properly  designed  optimal  rate  estimator  Kalman 
filter  is  effective  and  able  to  estimate  body  rate  from  a  single  star  sensor.  In  addition, 
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these  initial  results  prove  that  a  single  sensor  couple  with  a  proper  rate  estimator  design 
can  be  used  as  backup  or  even  primary  attitude  determination  process. 

For  NPSAT-1,  a  single  star  sensor  estimator  will  be  an  addition  to  the  control 
system.  This  approach  will  be  necessary,  especially  during  the  initial  launch  from  the 
Pegasus,  where  after  initial  launch  the  spacecraft  will  tumbling  at  some  pitch,  roll,  and 
yaw  rates.  The  magnetometer  and  the  magnetic  torquers  control,  but  would  require,  an 
additional  vector  which  the  star  sensor  can  provide 
C.  FUTURE  RESEARCH  AREAS 

The  future  of  small  satellite  design  and  development  is  an  interesting  area.  One 
approach  would  be  to  implement  a  MIMO  controller  using  the  state  space  approach.  This 
could  allow  for  more  sophisticated  techniques  and  possibly  lower  the  design 
requirements  of  the  sensors  and  actuators. 

A  second  approach  might  be  to  design  a  controller  that  would  eliminate  the 
momentum  wheel,  replacing  it  with  only  three-axis  magnetic  torquers.  This  advance 
requires  extensive  investigation  and  simulation.  However,  initial  testing  indicates  that 
this  may  be  an  adequate  design  approach  for  some  future  mission  with  less  stringent 
pointing  accuracy.  Additional  simulation  would  be  necessary  using  a  Monte  Carlo 
approach  for  different  as  and  /us  as  required  (a  Monte  Carlo  approach  is  considered 
suitable). 

Thirdly,  new  research  on  the  development  of  an  optimal  estimator  that  includes 
‘sensor  fusion.’  Sensor  fusion  would  combine  all  rates  from  every  sensor  and  merge  the 
data  into  a  more  accurate  estimator.  Due  to  the  failure  of  aging  satellites,  design  of  a 
sensor  fusion  type  of  rate  estimator  could  be  useful  to  replace  in-orbit,  failing  rate  gyros, 
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and  allow  for  updating  software  packages  that  perform  the  rate  estimation  for  gyros. 
Data  streams  would  emanate  from  several  different  types  of  onboard  sensor  devices  such 
as:  star  trackers,  horizon  sensors,  magnetometers,  and  sun  sensors 

A  fourth  area  of  research  would  involve  the  concept  of  reverse  time  smoothing. 
Typically,  the  Kalman  Filter  uses  only  past  and  present  observations,  and  is  therefore  a 
causal  filter.  This  is  ideal  for  real  time  systems  such  as  satellites.  However,  for 
improved  estimates,  the  additional  computing  power  of  modern  satellites  could  be  used  to 
post-process  old  data.  The  smoothed  past  estimates  could  then  be  used  in  the  Kalman 
Filter. 

When  considering  a  fixed  interval  smoother,  several  methods  have  been 
developed  to  post-process  data.  Three  of  the  most  common  are  1)  forward-backward 
smoother,  2)  two-point  boundary  value  approach,  and  3)  the  Rauch-Tung-Streibel 
smoother.  Additional  work  in  this  area  could  potentially  offer  some  accuracy 
improvements. 


71 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


72 


Appendix  Figure  1 :  Simulink  Block  Diagram. 
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Appendix  Figure  4:  H-Block  Integrator. 


Large  angles 


Appendix  Figure  5:  Rotation  Integrator. 


Appendix  Figure  6:  Rotation  Integrator  E-Block. 
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Appendix  Figure  7:  Rotation  Integrator  e2b_C_o. 
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Appendix  Figure  8:  w  n  body. 
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Appendix  Figure  9:  b_C_o  to  phi  theta  psi. 
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function  M=Dumping(u) 

%wl=u ( 1 ) ; 
w2=u (2 ) ; 
cl3=u  (3)  ; 

%c2  3=u  (4)  ; 

Kl=14 ; 

K2=0 . 1146; 

K3=l ; 

T=K2  * (-Kl*w2+cl3) ; 

M=K3*T* [0  1  0]  '  ; 

function  b_C_o=e2b_C_o (e) 
e=reshape (e, 4 , 1 ) ; 

C=zeros (3,3) ; 

C(l,l)=e(l)*e(l)-e(2)*e(2)-e(3)*e(3)+e(4)*e(4); 

C  (2, 1)  =2*  (e(l)*e(2)-e(3)*e(4)  )  ; 

C  (3, 1)  =2*  (e(l)*e(3)+e(2)*e(4)  )  ; 

C  (1, 2)  =2*  (e(l)*e(2)+e(3)*e(4)  )  ; 

C  (2, 2)=e  (2)  *e  (2)  -e  (1)  *e  (1)  -e  (3)  *e  (3)  +e  (4)  *e  (4)  ; 
C(3,2)=2*(e(2)*e(3)-e(l)*e(4)  )  ; 

C  (1, 3)  =2*  (e(l)*e(3)-e(2)*e(4)  )  ; 

C(2,3)=2*  (e(2)*e(3)+e(l)  *e(4)  )  ; 

C(3, 3)  =e  (3)  *e  (3)  -e  (1)  *e  (1)  -e  (2)  *e  (2)  +e  (4)  *e  (4)  ; 

Error=0 ; 

tol=le-6; 

if  (abs (norm (C ( 1 , : ) ) -1 ) >tol ) 

Error=l ; 

elseif  (abs (norm (C (2 , : ) ) -1 ) >tol) 

Error=l ; 

elseif  (abs (norm (C (3, : ) ) -1 ) >tol) 

Error=l ; 

elseif  (abs (norm (C (:, 1 )) -1 ) >tol ) 

Error=l ; 

elseif  (abs (norm (C ( : , 2 ) ) -1 ) >tol) 

Error=l ; 

elseif  (abs (norm (C ( : , 3) ) -1 ) >tol) 

Error=l ; 
end; 

if  (Error==l) 

disp(' Error  on  b  C  o! ! (Reduce  time  step  or  tighten 
tolerance . ' ) ; 

end; 

b_C_o=reshape (C,  9,1)  ; 

function  M=GravityGrad (u) 
b_C_o=reshape (u,  3,3)  ; 

C=b_C_o ( : , 3 ) ; 

I'  [22.63  0  0;  0  24.67  0;0  0  11]; 

%I= [2246.87  0  0;  0  3202.94  0;  0  0  969.121]; 
wo= .0010949; 

M=3*woA2*cross (C, I*C) ; 

function  B=MagneticField (U) 
b_C_o=U (1 : 9) ; 
t=U (10)  ; 
alpha0=0 ; 


79 


u0=0;i=75*pi/180; 

e=ll*pi/180; 

K=7 . 943el5; 
wO=. 0010949; 
we=2*pi/ (24*3600) ; 

C=reshape (b_C_o, 3,3)  ; 

alpha=alphaO+wO*t; 

u=u0+we*t; 

R=(6378. 4+550) *1000; 
r_o= [ 0  0  -1] ' ; 

Ca=cos (alpha) ; Sa=sin (alpha) ; 

Cu=cos (u) ; Su=sin  (u) ; 

Ce=cos (e) ; Se=sin  (e) ; 

Ci=cos ( i ) ; Si=sin  ( i )  ; 

Cx=-Ca* (Ce*Si-Se*Ci*Cu) +Sa*Se*Su; 

Cy=  Ce*Ci+Se*Si*Cu; 

Cz=  Se* (Ca*Su-Sa*Cu*Ci) +Ce*Sa*Si; 

M_o= [Cx, Cy, Cz ] ' ; 

B_o= (K/RA3) * (3* (M_o ' *r_o) *r_o-M_o) ; 

B_b=C*B_o; 

function  M=ReactionJets (u) 
wl=u  ( 1 ) ; 

%w2=u (2 ) ; 

%cl3=u  (3)  ; 
c2  3=u  (4)  ; 

Kl=  800; 

K2=  0.0573; 

K3=  0.268; 

K4  =  l ; 

Mx=-K2* (800*wl+c23) ; 

Mz=-K3*Mx; 

M= [Mx, 0 ,  Mz  ]  '  ; 

function  edot=we2edot (u) 

w=u (1:3) ; 

e=u (4:7) ; 

w=reshape (w, 3,1) ; 

e=reshape (e, 4 , 1 )  ; 

if  (abs (norm (e) -1 ) >le-6) 

disp('Sum  of  squares  of  e' '  not  summing  to  1 
end; 


M=  [  0 

w  (3) 

-w  (2) 

w  ( 1 )  ; 

-w  (3) 

0 

w  (1) 

w  ( 2  )  ; 

w  (2) 

-w  (1) 

0 

w  ( 3 )  ; 

-w  (1) 

-w  (2) 

-w  (3) 

0]  ; 

edot=0 . 5*M*e; 

function  wo=WoInBody (b  C  o) 
C=reshape (b_C_o, 3,3) ; 
wo= .0010949; 
wo=-wo*C ( : , 2 ) ; 

function  wo=Wo!nBody (b_C_o) 
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C=reshape (b_C_o, 3,3) ; 
wo= .0010949; 
wo=-wo*C (  :  ,  2  )  ; 


%  Sat  att  rate  est.. 

%  xl  is  roll,x2  is  roll  rate,  x3  is  yaw,  x4  is  yaw  rate, 
xp  is  pitch 

%  orbit  by  discrete.  This  becomes  the  pitch  observable 
xo  ( 3 , i ) 

Ao=  [  0  100;  0000;  0001;  0000] 

Bo= [  0  0 ; 1  0 ; 0  0;0  1] 

Co= [1  0  0  0 ]  ; 

%N=-1 ; Tf=2  94  00 ; dt=2 ; kmax=Tf / dt  +1 ; 

N=-l ;Tf=38900; dt=2 ; kmax=Tf / dt  +1 ; 

[Phio,  Delo] =c2d (Ao, Bo, dt) 

uo=zeros (2,  kmax) ;ho=zeros (1, kmax) ;xo=zeros (4, kmax) ; 
yo=zeros ( 1 , kmax) ; time=zeros ( 1 , kmax) ; xcart=zeros (2 , kmax) ; 
xo(;,l)=[7439*0.62  0  0  0.00130]'; 
for  (i=l:kmax-l) 

uo(:,i)=[xo(l,i) *xo(4,i) A2- 

(32 .16/5280) *4000A2/ (xo (1, i) ) A2; 

- (2*xo (2, i) *xo (4, i) ) /  (xo  (1, i) ) ]  ; 

xo(:,i+l)  =  Phio*xo(:,i)  +Delo*uo(:,i); 
time(i+l)=  time(i)  +  dt; 
end; 

for  ( i=l : kmax) 

xcart (1, i) =xo (l,i)*sin(xo(3,i)  )  ; 
xcart (2, i) =xo (l,i)*cos(xo(3,i)  )  ; 
end 

w=0 .00013; 

A=[0  1  0  0;0  0  0  -w;0  0  0  1;0  w  0  0]  ; 

B=  [  0  0 ;  1  0 ; 0  0 ; 0  1 ]  ; 

Bp=  [  0  1]  '  ; 
dt=2  ; 

pz= [0.78  0.79  0.77  0.76]'; 

[Phi,  Del]  =c2d  (A,  B,  dt) 
k=place ( Phi , Del , pz ) 
ppp=eig ( Phi ) 
pppp=eig (Phi-Del*k) 
pause 

[Phi,  Del]  =c2d  (A,  B,  dt) 

Ap=  [  0  1;  0  0] 

[Phip,  Delp] =c2d (Ap, Bp, dt) 
p= [ . 85  . 856]  ' 
kd=place (Phip, Delp, p) 
pp=eig (Phip+Delp*kd) 

%pause 

u=zeros (2, kmax) ; 
up=zeros ( 1 , kmax) ; 
x=zeros ( 4 , kmax) ; 
xp=zeros (2, kmax) ; 
y=zeros ( 1 , kmax) ; 
xf=zeros (4, kmax) ; 
time=zeros ( 1 , kmax) ; 
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zkkml  =  zeros  (2,1)  ; 

%  Random  number  generator 
v=randn ( 1 , v) ; 
zkkmlp=0 ; 

P=50*eye (4) ; 

Q=[(dtA4)/4  (dtA3)/2  0  0; (dtA3) /2  dtA2  00;  ... 

0  0  (dtA4)/4  ( dt A 3 ) / 2 ; 0  0  (dtA3) /2  dt A2 ] *0 . 0001 ; 

R=  [ .0002  0 ; 0  .0002]; 

H= [ 1  0  0  0 ; 0  0  1  0  ]  ; 

Pp=50*eye (2 ) ; 

Qp=[(dtA4)/4  (dtA3) /2; (dtA3) /2  dtA2]*0.0001 
Rp= [. 00002 ] ; 

Hp= [ 1  0  ] ; 

x(:,l)=[0.001;  0.002  ;0.003;  0.004]; 

xkk=zeros (4, kmax) ; 

xkkml=zeros (4, kmax) ; 

xkkml (:,1)  =  [0.0;0.0;0;  0.0013]  ; 

xp (2, 1) =0 . 0013; 

xkkp=zeros (2 , kmax) ; 

xkkmlp=zeros (2, kmax) ; 

xkkmlp ( : , 1 )  =  [  0 ; 0 ]  ; 

zkkml  ( : ,  1)  =*j[0;0]  ; 

z=  [x  (1,  1)  ;x  (3,  1)  ]  ; 

for  (i=l : kmax-1) 

%  plant 

u  (1,  i)  =  - k  ( 1 ,  : )  * x  ( : ,  i )  ; 

u (2, i) =  -k(2,  : ) *  x ( : , i ) ; 

x(:,i+l)  =  Phi*x(:,i)  +  Del*u(:,i)  ; 

time(i+l)=  time(i)  +  dt; 

%  Random  number  generator 
v=randn ( 1 , v) ; 

y (1,  i  +  1) =0 . 00000 l*randn (1) ; 

%  Kalman  Filter 

G=P*H ' *inv (H*P*H ' +R) ; 

Pk= (eye (4) -G*H) *P; 

P=Phi*Pk*Phi ' +Q; 

xkk ( : , i) =xkkml ( : , i) +G* (z-zkkml) ; 
xkkml ( : , i+1) =Phi*xkk ( : ,  i )  ; 
zkkml = [xkkml ( 1 , i+1 ) ; xkkml (3,i+l)]; 
z= [x (1, i  +  1) ;x (3, i  +  1)  ]  ; 

%  from  Kepler  H  =  mrlA2*wl  =  mr2 A2*w2 . . . . r2  =  xo(2,i) 
(rlA2*wl/w2) A (0.5) 

%r2=( (7439*0.62) A2*0. 0013) A0. 5; 

Gp=Pp*Hp ' *inv (Hp*Pp*Hp ' +Rp) ; 

Pkp= (eye (2 ) -Gp*Hp) *Pp; 

Pp=Phip*Pkp*Phip ' +Qp; 
zp= [xo (3, i) +y (1, i) ] ; 

xkkp ( : , i) =xkkmlp ( : , i) +Gp* (zp-zkkmlp) ; 

up ( i ) =  - (2*xo (2, i) *xkkp (2, i) ) / (xo (1,  i) )  ; 
xp(:,i+l)  =  Phip*xp(:,i)  +Delp*up(i); 
xkkmlp ( : , i+1 ) =Phip*xkkp ( : , i ) +  Delp*up ( i ) ; 
zkkmlp= [xkkmlp ( 1 ,  i  +  1 ) ]  ; 

end 
cl  f 

figure  ( 1 ) 

plot (time (1, : ) , y (1, : ) ) 
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%pause 
%c  1  f 

figure  (2 ) 

%subplot (221) , 

plot (x (1, : ) ,  x  (2,  : )  ) 

title ('roll  vs  roll  rate  ') 

xlabel ( ' xl ' ) ,  ylabel ( ' x2 ' ) , grid 

%subplot (222 ) , 

figure ( 3 ) 

plot (time (1, : ) , xkk (3, : ) , ' o ' , time (l,:),x(3,:)) 
title ('yaw  est  vs  time') 

%subplot (223) , 
figure  ( 4 ) 

plot (x (3, : ) ,  x  (4,  : )  ) 
title ( '  yaw  vs  yaw  rate') 
xlabel ( ' x3 ' ) ,  ylabel ( ' x4 ' ) , grid 
%subplot (224  ) , 
figure ( 5 ) 

plot (time (1, : ) , xkk (4, : ) ) 

title ('yaw  rate  est  vs  time  ') 

xlabel ( ' xkk3 ' ) ,  ylabel ( ' xkk4 ' ) , grid 

%pause 

%clf  ; 

figure (6) 

%subplot (221) , 

plot (time ( 1 , : ) , xkk ( 1 , : ) , time ( 1 ,:), x ( 1 ,:),'.')  ; 
xlabel ( ' time ' ) , ylabel ( ' xkkl  &  xl '  )  ; 
title ( [ ' roll '  ] )  ; 

%subplot  (222 )  , 
figure  ( 7 ) 

plot (time ( 1 ,  : ) , xkkp (2 ,  : ) ,  '  x ' , time ( 1 ,  : ) , xp (2 ,  : )  ,  '  .  '  ) 
xlabel ( ' time ' ) , ylabel ( ' pitch  rate  xp2 ' ) ; 
title  (["]); 
grid; 

%subplot (223) , 
figure  ( 8 ) 

plot (time ( 1 ,  : ) , xkk (2 ,  : ) , time (1, : ) , x  (2  ,:),'.') 
xlabel ( ' time ' ) , ylabel ( ' roll  rate  xkk2 ' ) ; 

%subplot (224 ) , 
figure (9) 

plot (time ( 1 , : ) , xkkp (2 , : ) ) 

xlabel (' time '), ylabel (' pitch  rate  est  xkkp2 ' ) ; 
grid; 

%pause 

xxp= zeros (2,10) ; 

xxkkp=zeros (2, 10) ; 

ttime=zeros (1, 10)  ; 

for  ( i=l : 70 ) 

xxkk ( 2 , i ) =xkk ( 2  ,  i )  ; 

xx ( 2 , i ) =x ( 2 , i )  ; 

ttime ( 1 , i ) =time ( 1 , i ) ; 

end 

%clg 

figure (10) 

plot (time ( 1 , : ) , xp (2 , : ) , '  x ' , time ( 1 , : ) , xo ( 4 , : ) , ' o ' ) 
xlabel (' time '), ylabel (' pitch  rate  xo2  &  xp2 ' ) ; 
title  (["]); 
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%pause 
figure  (11) 

plot (xcart {!,:), xcart (2 , : ) ) 

title ( ' orbit  ' ) ; xlabel ( ' x ' ) ,  ylabel ( ' y ' ) , grid; 

%pause 
figure ( 12 ) 

plot (ttime ( 1 , : ) , xxkk (2 , : ) , '  x ' , ttime ( 1 , : ) , xx (2 , : ) , ' o ' ) 
xlabel ( ' time ' ) , ylabel ( ' roll  rate  x2 ' ) ; 
title ( [ ' ']); 

%pause 
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% MAT LAB  CODE 
%Primary  Code 

%  Code  written  by  Professor  Hal  Titus 
%  Modified  by  LCDR  John  Vitalich  01  JUN  02 

%  orbit  by  discrete.  This  becomes  the  pitch  observable  xo(3,i) 
%  from  Kepler  H  =  mrlA2*wl  =  mr2 A2*w2 . . . . r2  =  xo(2,i)  = 
(rlA2*wl/w2) A (0.5) 

%  r2=( (7439*0.62) A2*0. 0013) A0. 5; 


%  constants 
mu  =  3.986e5 
w=0. 00013; 
Tf=29400*2; 
Tf=50 ; 
dt=l ; 


%  Gravitational  coefficient  (km3/sec2) 
%  orbital  frequency  (rads/sec) 

%  seconds  in  one  12  hour  pass 

%  seconds  in  one  12  hour  pass 


kmax=Tf/dt  +1; 
trackerr=6*pi/ (3600*180) ; 


%  Kalman  Filter  Covariances  and  observation  matrices 
%  roll  yaw 
P=1 *eye ( 4 ) ; 
q=0 . 01 

Q=[(dtA4)/4  (dtA3)/2  0  0; (dtA3) /2  dtA2  00;  ... 

0  0  (dtA4)/4  ( dt A 3 ) / 2 ; 0  0  (dtA3) /2  dtA2]*q; 

R=  [trackerr A2  0;0  trackerrA2]; 

H=  [  1  0  0  0 ;  0  0  1  0  ]  ; 

%  pitch 
Pp=50*eye (2 ) ; 

Qp=[(dtA4)/4  (dtA3) /2; (dtA3) /2  dtA2]*q; 

Rp= [trackerr A2 ] ; 

Hp= [ 1  0  ]  ; 


%  initialize  matrices 

uo= zeros (2, kmax) ;ho= zeros (1, kmax) ;xo= zeros (4, kmax) ; 

yo=zeros ( 1 , kmax) ; time=zeros ( 1 , kmax) ; xcart=zeros (2 , kmax) ; 

u= zeros (2, kmax) ;up= zeros (1, kmax) ;x= zeros (4, kmax) ;xp= zeros (2, kmax) ; 

y= zeros ( 1 , kmax) ;xf=zeros (4, kmax) ; time= zeros ( 1 , kmax) ; zkkml= zeros (2,1) 

xkk=zeros (4, kmax) ;xkkml  =  zeros (4, kmax) ;  z=zeros (2,1)  ; 

xkkp=zeros (2 , kmax) ; xkkmlp=zeros (2  ,  kmax)  ; 


%  Define  the  Molniya  orbit 

Ao= [0  1  0  0 ; 0  0  0  0;0  0  0  1 ; 0  0  0  0 ]  ; 

Bo= [  0  0 ; 1  0 ; 0  0 ;  0  1  ]  ; 

Co= [100  0]; 

[Phio,  Delo] =c2d (Ao, Bo, dt) ; 


%  xol  is  r  (km) ,  xo2  is  r  dot,  xo3  is  theta  (rad) ,  xo4  is  theta  dot 
xo(:,l)  =  [7439  0  0  .0013] 

h=xo ( 1 ) A2 *xo  ( 4 )  %  angular  momentum 

for  (i=l:kmax-l) 

uo ( : , i) = [xo (1, i) *xo (4, i) A2-mu/ (xo (1, i) ) A2; 

- (2*xo (2, i) *xo (4, i) ) / (xo  (1, i) ) ] ; 
xo(:,i+l)  =  Phio*xo(:,i)  +Delo*uo(:,i); 
time(i+l)=  time(i)  +  dt; 
end; 

%  convert  orbit  to  rectangular  co-ordinates  for  plotting 
for  ( i=l : kmax) 
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xcart (1, i) =xo (1, i) *sin(xo(3,i) ) ; 
xcart (2,  i) =xo (l,i)*cos(xo(3,i) ) ; 
xc(l,i)=  6378*sin (xo (3, i) ) ; 
xc(2,i)=  637 8 *cos (xo  ( 3 , i ) ) ; 
end 

%  xl  is  roll,x2  is  roll  rate,  x3  is  yaw,  x4  is  yaw  rate,  xp  is  pitch 
A=[0  1  0  0;0  0  0  -w;0  0  0  1;0  w  0  0] ; 

B=  [  0  0 ;  1  0 ; 0  0 ; 0  1 ]  ; 

Ap= [0  1 ;  0  0 ] ; 

Bp=  [  0  1]  '  ; 

%  place  desired  poles  for  roll  and  yaw 
pz= [0.78  0.79  0.77  0.76] ' ; 

%pz=[ .4  .41  .42  .43]  '  ; 

%  state  transition  matrix 
[Phi,  Del]  =c2d  (A,  B,  dt)  ; 

%  gains  and  eigenvalues 
k=place ( Phi , Del , pz ) ; 

[klqr, s, elqr] =dlqr (Phi, Del, Q, R)  %  discrete  linear  quadratic 

regulator 

%k=klqr ; 

ppp=eig ( Phi ) ; 

pppp=eig (Phi-Del*k) ; 

%  place  desired  poles  for  pitch 
p= [ . 85  . 856] ' ; 

[Phip,  Delp] =c2d (Ap, Bp, dt) ; 

%  gains  and  eigenvalues 
kd=place ( Phip, Delp, p) ; 

[klqrp, s, elqr] =dlqr (Phip, Delp, Qp, Rp)  %  discrete  linear 

quadratic  regulator 

%kd=klqrp; 

pp=eig (Phip+Delp*kd)  ; 


%  Kalman  filter 
%  initial  conditions 

x ( : , 1) = [0 . 001;  0.002  ;0.003;  0.004];  %  x  is  truth  state 

xkkml ( : , 1) = [0 . 0; 0 . 0; 0; 0 . 0013] ;  %  z  is  truth  measurement 

% (includes  meas .  error) 

zkkmlp=0 ; 

xp (2, 1) =0 . 0013; 

xkkmlp ( ; , 1 ) = [ 0 ; 0 ] ; 

zkkml ( : , 1) = [ 0 ; 0 ] ; 

zp= .0013; 

%  run  filter 
for  (i=l:kmax-l) 

%  plant 

time(i+l)=  time(i)  +  dt; 


%  for  roll  and  yaw 
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G=P*H ' *inv (H*P*H ' +R)  ; 

Pk= (eye (4) -G*H)  *P; 

Update 

P=Phi*Pk*Phi ' +Q; 

Prediction 

xkk ( : , i) =xkkml ( : , i) +G* (z-zkkml)  ; 
u ( 1 , i ) =  k (1, : ) * (xkk ( : , i ) -x ( : , i ) ) ; 
u (2, i) =  k (2,  : ) * (xkk ( : , i ) -x  ( : ,  i ) )  ; 
x(:,i+l)  =  Phi*x(:,i)  +  Del*u(:,i)  ; 
xkkml ( : , i+1) =Phi*xkk ( : , i) +  Del*u(:,i)  ; 
estimate 

roller=12* (randn (1) - . 5) *pi/ (3600*180)  ; 
random  error 

yawer=12* (randn (1) -.5) *pi/ (3600*180)  ; 
random  error 

pitcher=12* (randn (1) -.5) *pi/ (3600*180) ; 
random  error 

z= [ 1+roller ; 2+yawer ] ; 
zkkml= [xkkml ( 1 , i+1 ) ; xkkml (3 , i+1 ) ] ; 
estimate 

%  for  pitch 

Gp=Pp*Hp ' *inv (Hp*Pp*Hp ' +Rp) ; 

Pkp= (eye (2 ) -Gp*Hp) *Pp; 

Update 

Pp=Phip*Pkp*Phip ' +Qp; 

Prediction 

xkkp ( : , i) =xkkmlp ( : , i) +Gp* (zp-zkkmlp)  ; 
update 

up ( i ) =  -  (2*xo (2, i) *xkkp (2, i) ) / (xo (1, i) )  ; 
xp(:,i+l)  =  Phip*xp(:,i)  +Delp*up(i); 
xkkmlp ( : , i+1 ) =Phip*xkkp ( : , i ) +  Delp*up ( i ) ; 
estimate 

zp= [xo ( 3 , i ) +pitcher ] ; 

Measurement 

zkkmlp= [xkkmlp ( 1 ,  i  +  1 ) ]  ; 

Measurement  estimate 
end  %  Kalman  loop 


%  Kalman  Gain 

%  Covariance 

%  Covariance 

%  State  estimate  update 
%  Control  Law 
%  Control  Law 
%  Update  State 
%  Predict  State 


%  +/-  6  arcseconds  of 
%  +/-  6  arcseconds  of 
%  +/-  6  arcseconds  of 


%  Update  Measurement 


%  Kalman  Gain 

%  Covariance 

%  Covariance 

%  State  estimate 

%  Control  Law 
%  Update  State 
%  Predict  State 

%  Update 

%  Update 


cl  f 

figure ( 1 ) 

plot (xcart (1,  : ) , xcart (2,  : ) ,0,0,  ' * ' ,xc  (1,  : ) , xc (2,  : ) ) 

AXIS ([-35000  35000  -60000  10000]) 

XLABEL ( ' km ' ) ;  YLABEL ( ' km ' ) 

figure (2 ) 

subplot (211) , plot (x ( 1 , :) ,x(2, :) ) 
title ('roll  vs  roll  rate  ') 
xlabel ( ' xl ' ) ,  ylabel ( ' x2 ' ) , grid 

subplot (212) , plot (time ( 1 , : ) , xkk ( 1 , : ) , ' . ' , time (1, :) ,x(l, :) ) 
title ('roll  &  roll  estimate  vs  time') 
xlabel ( ' time ' ) ,  ylabel ( ' x3  &  xkk3 ' ) , grid 
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figure ( 3 ) 

subplot (211) , plot (x ( 3 , :) , x ( 4 , :) ) 
title ( '  yaw  vs  yaw  rate') 
xlabel ( ' x3 ' ) ,  ylabel ( ' x4 ' ) , grid 

subplot (212) , plot (time (1 ,:), xkk (3 , time ( 1 , 
title ('yaw  &  yaw  estimate  vs  time') 
xlabel ( ' time ' ) ,  ylabel ( ' x3  &  xkk3 ' ) , grid 

figure  ( 4 ) 
subplot (211) 

plot (time ( 1 , : ) , xkkp ( 1 , : ) , time ( 1 , : ) , xp ( 1 , : ) , ' . ' ) 
title ([' pitch  &  pitch  est  vs  time']); 
xlabel ( ' time ' ) , ylabel ( ' pitch  rate  xp2 ' ) ; 
grid; 

subplot (212) 

plot (time ( 1 ,  : ) , xkkp (2 ,  : ) , time ( 1 ,  : ) , xp  (2 ,  : ) ,  '  .  ' ) 
title ([' pitch  rate  &  pitch  rate  est  vs  time']); 
xlabel ( ' time ' ) , ylabel ( ' pitch  rate  xp2 ' ) ; 
grid; 

%  Profile  Generator 

%function  xtrue=prof ile ( kmax) 

%xtrue 

%  Random  number  generator 
for  i=l : kmax 

y(l, i+1) =12* (randn (1) -.5) *pi/ (3600*180) ;  %  +  / 

random  error 

y (2, i+1) =12* (randn (1) -.5) *pi/ (3600*180) ;  %  +/ 

random  error 

xtrue ( : , i ) = [ 1 ; 0 ; 2 ; 0 ; 1+y ( 1 , i+1 ) ;  2+y(2,i+l)];  % 

end 


: )  ,  x  (3,  : )  ) 


6  arcseconds  of 
6  arcseconds  of 

profile  variable 
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